diff --git a/src/pylife/materiallaws/notch_approximation_law.py b/src/pylife/materiallaws/notch_approximation_law.py index fb7f4d29..941a04a5 100644 --- a/src/pylife/materiallaws/notch_approximation_law.py +++ b/src/pylife/materiallaws/notch_approximation_law.py @@ -108,170 +108,6 @@ class ExtendedNeuber(NotchApproximationLawBase): ''' - def _e_star(self, load): - """Compute the plastic corrected strain term e^{\ast} from the Neuber approximation - (eq. 2.5-43 in FKM nonlinear) - - e_star = L/K_p / E + (L/K_p / K')^(1/n') - """ - - corrected_load = load / self._K_p - return self._ramberg_osgood_relation.strain(corrected_load) - - def _d_e_star(self, load): - """Compute the first derivative of self._e_star(load) - - e_star = L/K_p / E + (L/K_p / K')^(1/n') - - de_star(L)/dL = d/dL[ L/K_p / E + (L/K_p / K')^(1/n') ] - = 1/(K_p * E) + tangential_compliance(L/K_p) / K_p - """ - return 1/(self.K_p * self.E) \ - + self._ramberg_osgood_relation.tangential_compliance(load/self.K_p) / self.K_p - - def _neuber_strain(self, stress, load): - """Compute the additional strain term from the Neuber approximation - (2nd summand in eq. 2.5-45 in FKM nonlinear) - - (L/sigma * K_p * e_star) - """ - - e_star = self._e_star(load) - - # bad conditioned problem for stress approximately 0 (divide by 0), use factor 1 instead - # convert data from int to float - if not isinstance(load, float): - load = load.astype(float) - # factor = load / stress, avoid division by 0 - factor = np.divide(load, stress, out=np.ones_like(load), where=stress!=0) - - return factor * self._K_p * e_star - - def _stress_implicit(self, stress, load): - """Compute the implicit function of the stress, f(sigma), - defined in eq.2.5-45 of FKM nonlinear - - f(sigma) = sigma/E + (sigma/K')^(1/n') - (L/sigma * K_p * e_star) - """ - - return self._ramberg_osgood_relation.strain(stress) - self._neuber_strain(stress, load) - - def _d_stress_implicit(self, stress, load): - """Compute the first derivative of self._stress_implicit - - df/dsigma - """ - - e_star = self._e_star(load) - return self._ramberg_osgood_relation.tangential_compliance(stress) \ - - load * self._K_p * e_star \ - * -np.power(stress, -2, out=np.ones_like(stress), where=stress!=0) - - def _delta_e_star(self, delta_load): - """Compute the plastic corrected strain term e^{\ast} from the Neuber approximation - (eq. 2.5-43 in FKM nonlinear), for secondary branches in the stress-strain diagram - """ - - corrected_load = delta_load / self._K_p - return self._ramberg_osgood_relation.delta_strain(corrected_load) - - def _d_delta_e_star(self, delta_load): - """Compute the first derivative of self._delta_e_star(load) - - delta_e_star = ΔL/K_p / E + 2*(ΔL/K_p / (2*K'))^(1/n') - = ΔL/K_p / E + 2*(ΔL/(2*K_p) / K')^(1/n') - - d_delta_e_star(ΔL)/dΔL = d/dΔL[ ΔL/K_p / E + 2*(ΔL/(2*K_p) / K')^(1/n') ] - = 1/(K_p * E) + 2*tangential_compliance(ΔL/(2*K_p)) / (2*K_p) - = 1/(K_p * E) + tangential_compliance(ΔL/(2*K_p)) / K_p - """ - return 1/(self.K_p * self.E) \ - + self._ramberg_osgood_relation.tangential_compliance(delta_load/(2*self.K_p)) / self.K_p - - def _neuber_strain_secondary(self, delta_stress, delta_load): - """Compute the additional strain term from the Neuber approximation (2nd summand in eq. 2.5-45 in FKM nonlinear)""" - - delta_e_star = self._delta_e_star(delta_load) - - # bad conditioned problem for delta_stress approximately 0 (divide by 0), use factor 1 instead - # convert data from int to float - if not isinstance(delta_load, float): - delta_load = delta_load.astype(float) - # factor = load / stress, avoid division by 0 - factor = np.divide(delta_load, delta_stress, out=np.ones_like(delta_load), where=delta_stress!=0) - - return factor * self._K_p * delta_e_star - - def _stress_secondary_implicit(self, delta_stress, delta_load): - """Compute the implicit function of the stress, f(sigma), defined in eq.2.5-46 of FKM nonlinear""" - - return self._ramberg_osgood_relation.delta_strain(delta_stress) - self._neuber_strain_secondary(delta_stress, delta_load) - - def _d_stress_secondary_implicit(self, delta_stress, delta_load): - """Compute the first derivative of self._stress_secondary_implicit - Note, the derivative of `self._ramberg_osgood_relation.delta_strain` is: - d/dΔsigma delta_strain(Δsigma) = d/dΔsigma 2*strain(Δsigma/2) - = 2*d/dΔsigma strain(Δsigma/2) = 2 * 1/2 * tangential_compliance(Δsigma/2) - = self._ramberg_osgood_relation.tangential_compliance(delta_stress/2) - """ - - delta_e_star = self._delta_e_star(delta_load) - - return self._ramberg_osgood_relation.tangential_compliance(delta_stress/2) \ - - delta_load * self._K_p * delta_e_star \ - * -np.power(delta_stress, -2, out=np.ones_like(delta_stress), where=delta_stress!=0) - - def _load_implicit(self, load, stress): - """Compute the implicit function of the stress, f(sigma), - as a function of the load, - defined in eq.2.5-45 of FKM nonlinear. - This is needed to apply the notch approximation law "backwards", i.e., - to get from stress back to load. This is required for the FKM nonlinear roughness & surface layer. - - f(L) = sigma/E + (sigma/K')^(1/n') - (L/sigma * K_p * e_star(L)) - """ - - return self._stress_implicit(stress, load) - - def _d_load_implicit(self, load, stress): - """Compute the first derivative of self._load_implicit - - f(L) = sigma/E + (sigma/K')^(1/n') - (L/sigma * K_p * e_star(L)) - - df/dL = d/dL [ -(L/sigma * K_p * e_star(L))] - = -1/sigma * K_p * e_star(L) - L/sigma * K_p * de_star/dL - - """ - - return -1/stress * self.K_p * self._e_star(load) \ - - load/stress * self.K_p * self._d_e_star(load) - - def _load_secondary_implicit(self, delta_load, delta_stress): - """Compute the implicit function of the stress, f(Δsigma), - as a function of the load, - defined in eq.2.5-46 of FKM nonlinear. - This is needed to apply the notch approximation law "backwards", i.e., - to get from stress back to load. This is required for the FKM nonlinear roughness & surface layer. - - f(ΔL) = Δsigma/E + 2*(Δsigma/(2*K'))^(1/n') - (ΔL/Δsigma * K_p * Δe_star(ΔL)) - - """ - - return self._stress_secondary_implicit(delta_stress, delta_load) - - def _d_load_secondary_implicit(self, delta_load, delta_stress): - """Compute the first derivative of self._load_secondary_implicit - - f(ΔL) = Δsigma/E + 2*(Δsigma/(2*K'))^(1/n') - (ΔL/Δsigma * K_p * Δe_star(ΔL)) - - df/dΔL = d/dΔL [ -(ΔL/Δsigma * K_p * Δe_star(ΔL))] - = -1/Δsigma * K_p * Δe_star(ΔL) - ΔL/Δsigma * K_p * dΔe_star/dΔL - - """ - - return -1/delta_stress * self.K_p * self._delta_e_star(delta_load) \ - - delta_load/delta_stress * self.K_p * self._d_delta_e_star(delta_load) - def stress(self, load, *, rtol=1e-4, tol=1e-4): '''Calculate the stress of the primary path in the stress-strain diagram at a given elastic-plastic stress (load), from a FE computation. @@ -425,15 +261,179 @@ def load_secondary_branch(self, delta_stress, rtol=1e-4, tol=1e-4): ''' - # self._stress_implicit(stress) = 0 - # f(sigma) = sigma/E + (sigma/K')^(1/n') - (L/sigma * K_p * e_star) = 0 - # => sigma/E + (sigma/K')^(1/n') = (L/sigma * K_p * e_star) - # => (sigma/E + (sigma/K')^(1/n')) / K_p * sigma = L * e_star(L) - # <=> self._ramberg_osgood_relation.strain(stress) / self._K_p * stress = L * e_star(L) + # self._stress_implicit(stress) = 0 + # f(sigma) = sigma/E + (sigma/K')^(1/n') - (L/sigma * K_p * e_star) = 0 + # => sigma/E + (sigma/K')^(1/n') = (L/sigma * K_p * e_star) + # => (sigma/E + (sigma/K')^(1/n')) / K_p * sigma = L * e_star(L) + # <=> self._ramberg_osgood_relation.strain(stress) / self._K_p * stress = L * e_star(L) + + delta_load = optimize.newton(func=self._load_secondary_implicit, x0=delta_stress, fprime=self._d_load_secondary_implicit, \ + args=([delta_stress]), rtol=rtol, tol=tol, maxiter=20) + return delta_load + + def _e_star(self, load): + """Compute the plastic corrected strain term e^{\ast} from the Neuber approximation + (eq. 2.5-43 in FKM nonlinear) + + e_star = L/K_p / E + (L/K_p / K')^(1/n') + """ + + corrected_load = load / self._K_p + return self._ramberg_osgood_relation.strain(corrected_load) + + def _d_e_star(self, load): + """Compute the first derivative of self._e_star(load) + + e_star = L/K_p / E + (L/K_p / K')^(1/n') + + de_star(L)/dL = d/dL[ L/K_p / E + (L/K_p / K')^(1/n') ] + = 1/(K_p * E) + tangential_compliance(L/K_p) / K_p + """ + return 1/(self.K_p * self.E) \ + + self._ramberg_osgood_relation.tangential_compliance(load/self.K_p) / self.K_p + + def _neuber_strain(self, stress, load): + """Compute the additional strain term from the Neuber approximation + (2nd summand in eq. 2.5-45 in FKM nonlinear) + + (L/sigma * K_p * e_star) + """ + + e_star = self._e_star(load) + + # bad conditioned problem for stress approximately 0 (divide by 0), use factor 1 instead + # convert data from int to float + if not isinstance(load, float): + load = load.astype(float) + # factor = load / stress, avoid division by 0 + factor = np.divide(load, stress, out=np.ones_like(load), where=stress!=0) + + return factor * self._K_p * e_star + + def _stress_implicit(self, stress, load): + """Compute the implicit function of the stress, f(sigma), + defined in eq.2.5-45 of FKM nonlinear + + f(sigma) = sigma/E + (sigma/K')^(1/n') - (L/sigma * K_p * e_star) + """ + + return self._ramberg_osgood_relation.strain(stress) - self._neuber_strain(stress, load) + + def _d_stress_implicit(self, stress, load): + """Compute the first derivative of self._stress_implicit + + df/dsigma + """ + + e_star = self._e_star(load) + return self._ramberg_osgood_relation.tangential_compliance(stress) \ + - load * self._K_p * e_star \ + * -np.power(stress, -2, out=np.ones_like(stress), where=stress!=0) + + def _delta_e_star(self, delta_load): + """Compute the plastic corrected strain term e^{\ast} from the Neuber approximation + (eq. 2.5-43 in FKM nonlinear), for secondary branches in the stress-strain diagram + """ + + corrected_load = delta_load / self._K_p + return self._ramberg_osgood_relation.delta_strain(corrected_load) + + def _d_delta_e_star(self, delta_load): + """Compute the first derivative of self._delta_e_star(load) + + delta_e_star = ΔL/K_p / E + 2*(ΔL/K_p / (2*K'))^(1/n') + = ΔL/K_p / E + 2*(ΔL/(2*K_p) / K')^(1/n') + + d_delta_e_star(ΔL)/dΔL = d/dΔL[ ΔL/K_p / E + 2*(ΔL/(2*K_p) / K')^(1/n') ] + = 1/(K_p * E) + 2*tangential_compliance(ΔL/(2*K_p)) / (2*K_p) + = 1/(K_p * E) + tangential_compliance(ΔL/(2*K_p)) / K_p + """ + return 1/(self.K_p * self.E) \ + + self._ramberg_osgood_relation.tangential_compliance(delta_load/(2*self.K_p)) / self.K_p + + def _neuber_strain_secondary(self, delta_stress, delta_load): + """Compute the additional strain term from the Neuber approximation (2nd summand in eq. 2.5-45 in FKM nonlinear)""" + + delta_e_star = self._delta_e_star(delta_load) + + # bad conditioned problem for delta_stress approximately 0 (divide by 0), use factor 1 instead + # convert data from int to float + if not isinstance(delta_load, float): + delta_load = delta_load.astype(float) + # factor = load / stress, avoid division by 0 + factor = np.divide(delta_load, delta_stress, out=np.ones_like(delta_load), where=delta_stress!=0) + + return factor * self._K_p * delta_e_star + + def _stress_secondary_implicit(self, delta_stress, delta_load): + """Compute the implicit function of the stress, f(sigma), defined in eq.2.5-46 of FKM nonlinear""" + + return self._ramberg_osgood_relation.delta_strain(delta_stress) - self._neuber_strain_secondary(delta_stress, delta_load) + + def _d_stress_secondary_implicit(self, delta_stress, delta_load): + """Compute the first derivative of self._stress_secondary_implicit + Note, the derivative of `self._ramberg_osgood_relation.delta_strain` is: + d/dΔsigma delta_strain(Δsigma) = d/dΔsigma 2*strain(Δsigma/2) + = 2*d/dΔsigma strain(Δsigma/2) = 2 * 1/2 * tangential_compliance(Δsigma/2) + = self._ramberg_osgood_relation.tangential_compliance(delta_stress/2) + """ + + delta_e_star = self._delta_e_star(delta_load) + + return self._ramberg_osgood_relation.tangential_compliance(delta_stress/2) \ + - delta_load * self._K_p * delta_e_star \ + * -np.power(delta_stress, -2, out=np.ones_like(delta_stress), where=delta_stress!=0) + + def _load_implicit(self, load, stress): + """Compute the implicit function of the stress, f(sigma), + as a function of the load, + defined in eq.2.5-45 of FKM nonlinear. + This is needed to apply the notch approximation law "backwards", i.e., + to get from stress back to load. This is required for the FKM nonlinear roughness & surface layer. + + f(L) = sigma/E + (sigma/K')^(1/n') - (L/sigma * K_p * e_star(L)) + """ + + return self._stress_implicit(stress, load) + + def _d_load_implicit(self, load, stress): + """Compute the first derivative of self._load_implicit + + f(L) = sigma/E + (sigma/K')^(1/n') - (L/sigma * K_p * e_star(L)) + + df/dL = d/dL [ -(L/sigma * K_p * e_star(L))] + = -1/sigma * K_p * e_star(L) - L/sigma * K_p * de_star/dL + + """ + + return -1/stress * self.K_p * self._e_star(load) \ + - load/stress * self.K_p * self._d_e_star(load) + + def _load_secondary_implicit(self, delta_load, delta_stress): + """Compute the implicit function of the stress, f(Δsigma), + as a function of the load, + defined in eq.2.5-46 of FKM nonlinear. + This is needed to apply the notch approximation law "backwards", i.e., + to get from stress back to load. This is required for the FKM nonlinear roughness & surface layer. - delta_load = optimize.newton(func=self._load_secondary_implicit, x0=delta_stress, fprime=self._d_load_secondary_implicit, \ - args=([delta_stress]), rtol=rtol, tol=tol, maxiter=20) - return delta_load + f(ΔL) = Δsigma/E + 2*(Δsigma/(2*K'))^(1/n') - (ΔL/Δsigma * K_p * Δe_star(ΔL)) + + """ + + return self._stress_secondary_implicit(delta_stress, delta_load) + + def _d_load_secondary_implicit(self, delta_load, delta_stress): + """Compute the first derivative of self._load_secondary_implicit + + f(ΔL) = Δsigma/E + 2*(Δsigma/(2*K'))^(1/n') - (ΔL/Δsigma * K_p * Δe_star(ΔL)) + + df/dΔL = d/dΔL [ -(ΔL/Δsigma * K_p * Δe_star(ΔL))] + = -1/Δsigma * K_p * Δe_star(ΔL) - ΔL/Δsigma * K_p * dΔe_star/dΔL + + """ + + return -1/delta_stress * self.K_p * self._delta_e_star(delta_load) \ + - delta_load/delta_stress * self.K_p * self._d_delta_e_star(delta_load) class Binned: """Binning for notch approximation laws, as described in FKM nonlinear 2.5.8.2, p.55. @@ -459,109 +459,6 @@ def __init__(self, notch_approximation_law, maximum_absolute_load, number_of_bin self._number_of_bins = number_of_bins self._create_bins() - def _create_bins(self): - """Initialize the lookup tables by precomputing the notch approximation law values. - """ - # for multiple assessment points at once use a DataFrame with MultiIndex - if isinstance(self._maximum_absolute_load, pd.DataFrame): - assert self._maximum_absolute_load.index.name == "node_id" - - self._create_bins_multiple_assessment_points() - - # for a single assessment point use the standard data structure - else: - self._create_bins_single_assessment_point() - - - def _create_bins_single_assessment_point(self): - """Initialize the lookup tables by precomputing the notch approximation law values, - for the case of scalar variables, i.e., only a single assessment point.""" - - # create look-up table (lut) for the primary branch values, named PFAD in FKM nonlinear - self._lut_primary_branch = pd.DataFrame(0, - index=pd.Index(np.arange(1, self._number_of_bins+1), name="class_index"), - columns=["load", "strain", "stress"]) - - - self._lut_primary_branch.load \ - = self._lut_primary_branch.index/self._number_of_bins * self._maximum_absolute_load - - self._lut_primary_branch.stress \ - = self._notch_approximation_law.stress(self._lut_primary_branch.load) - - self._lut_primary_branch.strain \ - = self._notch_approximation_law.strain( - self._lut_primary_branch.stress, self._lut_primary_branch.load) - - # create look-up table (lut) for the secondary branch values, named AST in FKM nonlinear - # Note that this time, we used twice the number of entries with the same bin width. - self._lut_secondary_branch = pd.DataFrame(0, - index=pd.Index(np.arange(1, 2*self._number_of_bins+1), name="class_index"), - columns=["delta_load", "delta_strain", "delta_stress"]) - - self._lut_secondary_branch.delta_load \ - = self._lut_secondary_branch.index/self._number_of_bins * self._maximum_absolute_load - - self._lut_secondary_branch.delta_stress \ - = self._notch_approximation_law.stress_secondary_branch(self._lut_secondary_branch.delta_load) - - self._lut_secondary_branch.delta_strain \ - = self._notch_approximation_law.strain_secondary_branch( - self._lut_secondary_branch.delta_stress, self._lut_secondary_branch.delta_load) - - def _create_bins_multiple_assessment_points(self): - """Initialize the lookup tables by precomputing the notch approximation law values, - for the case of vector-valued variables caused by an assessment on multiple points at once.""" - - # name column "max_abs_load" - self._maximum_absolute_load.rename(columns={self._maximum_absolute_load.columns[0]: "max_abs_load"}, inplace=True) - - # create look-up table (lut) for the primary branch values, named PFAD in FKM nonlinear - index = pd.MultiIndex.from_product([np.arange(1, self._number_of_bins+1), self._maximum_absolute_load.index], - names = ["class_index", "node_id"]) - - self._lut_primary_branch = pd.DataFrame(0, index=index, columns=["load", "strain", "stress"]) - - # create cartesian product of class index and max load - a = pd.DataFrame({"class_index": np.arange(1, self._number_of_bins+1)}) - class_index_with_max_load = a.merge(self._maximum_absolute_load, how='cross') - - # calculate load of each bin - load = class_index_with_max_load.class_index.astype(float)/self._number_of_bins * class_index_with_max_load.max_abs_load - load.index = index - self._lut_primary_branch.load = load - - self._lut_primary_branch.stress \ - = self._notch_approximation_law.stress(self._lut_primary_branch.load) - - self._lut_primary_branch.strain \ - = self._notch_approximation_law.strain( - self._lut_primary_branch.stress, self._lut_primary_branch.load) - - # ---------------- - # create look-up table (lut) for the secondary branch values, named AST in FKM nonlinear - # Note that this time, we used twice the number of entries with the same bin width. - index = pd.MultiIndex.from_product([np.arange(1, 2*self._number_of_bins+1), self._maximum_absolute_load.index], - names = ["class_index", "node_id"]) - - self._lut_secondary_branch = pd.DataFrame(0, index=index, columns=["delta_load", "delta_strain", "delta_stress"]) - - # create cartesian product of class index and max load - a = pd.DataFrame({"class_index": np.arange(1, 2*self._number_of_bins+1)}) - class_index_with_max_load = a.merge(self._maximum_absolute_load, how='cross') - - # calculate load of each bin - delta_load = class_index_with_max_load.class_index.astype(float)/self._number_of_bins * class_index_with_max_load.max_abs_load - delta_load.index = index - self._lut_secondary_branch.delta_load = delta_load - - self._lut_secondary_branch.delta_stress \ - = self._notch_approximation_law.stress_secondary_branch(self._lut_secondary_branch.delta_load) - - self._lut_secondary_branch.delta_strain \ - = self._notch_approximation_law.strain_secondary_branch( - self._lut_secondary_branch.delta_stress, self._lut_secondary_branch.delta_load) - @property def ramberg_osgood_relation(self): '''The ramberg osgood relation object @@ -722,8 +619,7 @@ def strain(self, stress, load): raise ValueError(f"Binned class is initialized with a maximum absolute load of {self._maximum_absolute_load}, "\ f" but a higher absolute load value of |{load}| is requested (in strain()).") - return sign * self._lut_primary_branch.iloc[index+1].strain # "+1", because the next higher class is used - + return sign * self._lut_primary_branch.iloc[index+1].strain # "+1", because the next higher class is used def stress_secondary_branch(self, delta_load, rtol=1e-5, tol=1e-6): '''Get the stress on secondary branches in the stress-strain diagram at a given load @@ -879,4 +775,106 @@ def strain_secondary_branch(self, delta_stress, delta_load): f" but a higher absolute delta_load value of |{delta_load}| is requested (in strain_secondary_branch()).") return sign * self._lut_secondary_branch.iloc[index+1].delta_strain # "-1", transform to zero-based indices - \ No newline at end of file + + def _create_bins(self): + """Initialize the lookup tables by precomputing the notch approximation law values. + """ + # for multiple assessment points at once use a DataFrame with MultiIndex + if isinstance(self._maximum_absolute_load, pd.DataFrame): + assert self._maximum_absolute_load.index.name == "node_id" + + self._create_bins_multiple_assessment_points() + + # for a single assessment point use the standard data structure + else: + self._create_bins_single_assessment_point() + + def _create_bins_single_assessment_point(self): + """Initialize the lookup tables by precomputing the notch approximation law values, + for the case of scalar variables, i.e., only a single assessment point.""" + + # create look-up table (lut) for the primary branch values, named PFAD in FKM nonlinear + self._lut_primary_branch = pd.DataFrame(0, + index=pd.Index(np.arange(1, self._number_of_bins+1), name="class_index"), + columns=["load", "strain", "stress"]) + + + self._lut_primary_branch.load \ + = self._lut_primary_branch.index/self._number_of_bins * self._maximum_absolute_load + + self._lut_primary_branch.stress \ + = self._notch_approximation_law.stress(self._lut_primary_branch.load) + + self._lut_primary_branch.strain \ + = self._notch_approximation_law.strain( + self._lut_primary_branch.stress, self._lut_primary_branch.load) + + # create look-up table (lut) for the secondary branch values, named AST in FKM nonlinear + # Note that this time, we used twice the number of entries with the same bin width. + self._lut_secondary_branch = pd.DataFrame(0, + index=pd.Index(np.arange(1, 2*self._number_of_bins+1), name="class_index"), + columns=["delta_load", "delta_strain", "delta_stress"]) + + self._lut_secondary_branch.delta_load \ + = self._lut_secondary_branch.index/self._number_of_bins * self._maximum_absolute_load + + self._lut_secondary_branch.delta_stress \ + = self._notch_approximation_law.stress_secondary_branch(self._lut_secondary_branch.delta_load) + + self._lut_secondary_branch.delta_strain \ + = self._notch_approximation_law.strain_secondary_branch( + self._lut_secondary_branch.delta_stress, self._lut_secondary_branch.delta_load) + + def _create_bins_multiple_assessment_points(self): + """Initialize the lookup tables by precomputing the notch approximation law values, + for the case of vector-valued variables caused by an assessment on multiple points at once.""" + + # name column "max_abs_load" + self._maximum_absolute_load.rename(columns={self._maximum_absolute_load.columns[0]: "max_abs_load"}, inplace=True) + + # create look-up table (lut) for the primary branch values, named PFAD in FKM nonlinear + index = pd.MultiIndex.from_product([np.arange(1, self._number_of_bins+1), self._maximum_absolute_load.index], + names = ["class_index", "node_id"]) + + self._lut_primary_branch = pd.DataFrame(0, index=index, columns=["load", "strain", "stress"]) + + # create cartesian product of class index and max load + a = pd.DataFrame({"class_index": np.arange(1, self._number_of_bins+1)}) + class_index_with_max_load = a.merge(self._maximum_absolute_load, how='cross') + + # calculate load of each bin + load = class_index_with_max_load.class_index.astype(float)/self._number_of_bins * class_index_with_max_load.max_abs_load + load.index = index + self._lut_primary_branch.load = load + + self._lut_primary_branch.stress \ + = self._notch_approximation_law.stress(self._lut_primary_branch.load) + + self._lut_primary_branch.strain \ + = self._notch_approximation_law.strain( + self._lut_primary_branch.stress, self._lut_primary_branch.load) + + # ---------------- + # create look-up table (lut) for the secondary branch values, named AST in FKM nonlinear + # Note that this time, we used twice the number of entries with the same bin width. + index = pd.MultiIndex.from_product([np.arange(1, 2*self._number_of_bins+1), self._maximum_absolute_load.index], + names = ["class_index", "node_id"]) + + self._lut_secondary_branch = pd.DataFrame(0, index=index, columns=["delta_load", "delta_strain", "delta_stress"]) + + # create cartesian product of class index and max load + a = pd.DataFrame({"class_index": np.arange(1, 2*self._number_of_bins+1)}) + class_index_with_max_load = a.merge(self._maximum_absolute_load, how='cross') + + # calculate load of each bin + delta_load = class_index_with_max_load.class_index.astype(float)/self._number_of_bins * class_index_with_max_load.max_abs_load + delta_load.index = index + self._lut_secondary_branch.delta_load = delta_load + + self._lut_secondary_branch.delta_stress \ + = self._notch_approximation_law.stress_secondary_branch(self._lut_secondary_branch.delta_load) + + self._lut_secondary_branch.delta_strain \ + = self._notch_approximation_law.strain_secondary_branch( + self._lut_secondary_branch.delta_stress, self._lut_secondary_branch.delta_load) + \ No newline at end of file diff --git a/src/pylife/materiallaws/notch_approximation_law_seegerbeste.py b/src/pylife/materiallaws/notch_approximation_law_seegerbeste.py index 8d692af2..2279a85b 100644 --- a/src/pylife/materiallaws/notch_approximation_law_seegerbeste.py +++ b/src/pylife/materiallaws/notch_approximation_law_seegerbeste.py @@ -65,6 +65,201 @@ class SeegerBeste(pylife.materiallaws.notch_approximation_law.NotchApproximation used, which is a bit slower but arrives at a correct solution. ''' + def stress(self, load, rtol=1e-4, tol=1e-4): + '''Calculate the stress of the primary path in the stress-strain diagram at a given + elastic-plastic stress (load), from a FE computation. + This is done by solving for the root of f(sigma) in eq. 2.8-42 of FKM nonlinear. + + The secant method is used which does not rely on a derivative and has good numerical stability properties, + but is slower than Newton's method. The algorithm is implemented in scipy for multiple values at once. + The documentation states that this is faster for more than ~100 entries than a simple loop over the + individual values. + + We employ the scipy function on all items in the given array at once. + Usually, some of them fail and we recompute the value of the failed items afterwards. + Calling the Newton method on a scalar function somehow always converges, while calling + the Newton method with same initial conditions on the same values, but with multiple at once, fails sometimes. + + Parameters + ---------- + load : array-like float + The load + + Returns + ------- + stress : array-like float + The resulting stress + ''' + # initial value as given by correction document to FKM nonlinear + x0 = load * (1 - (1 - 1/self._K_p)/1000) + + # suppress the divergence warnings + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + + stress = optimize.newton(func=self._stress_implicit, x0=x0, \ + args=([load]), rtol=rtol, tol=tol, full_output=True, maxiter=50) + + # only for multiple points at once, if some points diverged + if hasattr(stress[1], "iloc") and sum(stress[1]) < len(stress[1]): + stress = self._stress_fix_not_converged_values(stress, load, x0, rtol, tol) + + return stress[0] + + def strain(self, stress, load): + '''Calculate the strain of the primary path in the stress-strain diagram at a given stress and load. + The formula is given by eq. 2.8-39 of FKM nonlinear. + load / stress * self._K_p * e_star + + Parameters + ---------- + stress : array-like float + The stress + load : array-like float + The load + + Returns + ------- + strain : array-like float + The resulting strain + ''' + + if not isinstance(stress, float): + stress = stress.astype(float) + + return self._ramberg_osgood_relation.strain(stress) + + def load(self, stress, rtol=1e-4, tol=1e-4): + '''Apply the notch-approximation law "backwards", i.e., compute the linear-elastic stress (called "load" or "L" in FKM nonlinear) + from the elastic-plastic stress as from the notch approximation. + This backward step is needed for the pfp FKM nonlinear surface layer & roughness. + + This method is the inverse operation of "stress", i.e., L = load(stress(L)) and S = stress(load(stress)). + + Note that this method is only implemented for the scalar case, as the FKM nonlinear surface layer & roughness + also only handles the scalar case with one assessment point at once, not with entire meshes. + + Parameters + ---------- + stress : array-like float + The elastic-plastic stress as computed by the notch approximation + + Returns + ------- + load : array-like float + The resulting load or linear elastic stress. + + ''' + + x0 = stress / (1 - (1 - 1/self._K_p)/1000) + + # suppress the divergence warnings + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + + load = optimize.newton(func=self._load_implicit, x0=x0, \ + args=([stress]), rtol=rtol, tol=tol, maxiter=50) + + return load + + def stress_secondary_branch(self, delta_load, rtol=1e-4, tol=1e-4): + '''Calculate the stress on secondary branches in the stress-strain diagram at a given + elastic-plastic stress (load), from a FE computation. + This is done by solving for the root of f(sigma) in eq. 2.8-43 of FKM nonlinear. + + Parameters + ---------- + delta_load : array-like float + The load increment of the hysteresis + + Returns + ------- + delta_stress : array-like float + The resulting stress increment within the hysteresis + + Todo + ---- + + In the future, we can evaluate the runtime performance and try a Newton method instead + of the currently used secant method to speed up the computation. + + .. code:: + + fprime=self._d_stress_secondary_implicit_numeric + + ''' + + # initial value as given by correction document to FKM nonlinear + x0 = delta_load * (1 - (1 - 1/self._K_p)/1000) + + # suppress the divergence warnings + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + + delta_stress = optimize.newton(func=self._stress_secondary_implicit, x0=x0, \ + args=([delta_load]), rtol=rtol, tol=tol, full_output=True, maxiter=50) + + # only for multiple points at once, if some points diverged + if hasattr(delta_stress[1], "iloc") and sum(delta_stress[1]) < len(delta_stress[1]): + delta_stress = self._stress_secondary_fix_not_converged_values(delta_stress, delta_load, x0, rtol, tol) + + return delta_stress[0] + + def strain_secondary_branch(self, delta_stress, delta_load): + '''Calculate the strain on secondary branches in the stress-strain diagram at a given stress and load. + The formula is given by eq. 2.8-43 of FKM nonlinear. + + Parameters + ---------- + delta_sigma : array-like float + The stress increment + delta_load : array-like float + The load increment + + Returns + ------- + strain : array-like float + The resulting strain + ''' + + if not isinstance(delta_stress, float): + delta_stress = delta_stress.astype(float) + + return self._ramberg_osgood_relation.delta_strain(delta_stress) + + def load_secondary_branch(self, delta_stress, rtol=1e-4, tol=1e-4): + '''Apply the notch-approximation law "backwards", i.e., compute the linear-elastic stress (called "load" or "L" in FKM nonlinear) + from the elastic-plastic stress as from the notch approximation. + This backward step is needed for the pfp FKM nonlinear surface layer & roughness. + + This method is the inverse operation of "stress", i.e., L = load(stress(L)) and S = stress(load(stress)). + + Note that this method is only implemented for the scalar case, as the FKM nonlinear surface layer & roughness + also only handles the scalar case with one assessment point at once, not with entire meshes. + + Parameters + ---------- + delta_stress : array-like float + The increment of the elastic-plastic stress as computed by the notch approximation + + Returns + ------- + delta_load : array-like float + The resulting load or linear elastic stress. + + ''' + + x0 = delta_stress / (1 - (1 - 1/self._K_p)/1000) + + # suppress the divergence warnings + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + + delta_load = optimize.newton(func=self._load_secondary_implicit, x0=x0, \ + args=([delta_stress]), rtol=rtol, tol=tol, maxiter=20) + + return delta_load + def _e_star(self, load): """Compute the plastic corrected strain term e^{\ast} from the Neuber approximation (eq. 2.5-43 in FKM nonlinear) @@ -222,7 +417,6 @@ def _d_stress_secondary_implicit_numeric(self, delta_stress, delta_load): h = 1e-4 return (self._stress_secondary_implicit(delta_stress+h, delta_load) - self._stress_secondary_implicit(delta_stress-h, delta_load)) / (2*h) - def _load_implicit(self, load, stress): """Compute the implicit function of the stress, f(sigma), as a function of the load, @@ -243,47 +437,6 @@ def _load_secondary_implicit(self, delta_load, delta_stress): return self._stress_secondary_implicit(delta_stress, delta_load) - def stress(self, load, rtol=1e-4, tol=1e-4): - '''Calculate the stress of the primary path in the stress-strain diagram at a given - elastic-plastic stress (load), from a FE computation. - This is done by solving for the root of f(sigma) in eq. 2.8-42 of FKM nonlinear. - - The secant method is used which does not rely on a derivative and has good numerical stability properties, - but is slower than Newton's method. The algorithm is implemented in scipy for multiple values at once. - The documentation states that this is faster for more than ~100 entries than a simple loop over the - individual values. - - We employ the scipy function on all items in the given array at once. - Usually, some of them fail and we recompute the value of the failed items afterwards. - Calling the Newton method on a scalar function somehow always converges, while calling - the Newton method with same initial conditions on the same values, but with multiple at once, fails sometimes. - - Parameters - ---------- - load : array-like float - The load - - Returns - ------- - stress : array-like float - The resulting stress - ''' - # initial value as given by correction document to FKM nonlinear - x0 = load * (1 - (1 - 1/self._K_p)/1000) - - # suppress the divergence warnings - with warnings.catch_warnings(): - warnings.simplefilter("ignore") - - stress = optimize.newton(func=self._stress_implicit, x0=x0, \ - args=([load]), rtol=rtol, tol=tol, full_output=True, maxiter=50) - - # only for multiple points at once, if some points diverged - if hasattr(stress[1], "iloc") and sum(stress[1]) < len(stress[1]): - stress = self._stress_fix_not_converged_values(stress, load, x0, rtol, tol) - - return stress[0] - def _stress_fix_not_converged_values(self, stress, load, x0, rtol, tol): '''For the values that did not converge in the previous vectorized call to optimize.newton, call optimize.newton again on the scalar value. This usually finds the correct solution.''' @@ -298,106 +451,7 @@ def _stress_fix_not_converged_values(self, stress, load, x0, rtol, tol): if result[1].converged: stress[0].iloc[index_diverged] = result[0] return stress - - def strain(self, stress, load): - '''Calculate the strain of the primary path in the stress-strain diagram at a given stress and load. - The formula is given by eq. 2.8-39 of FKM nonlinear. - load / stress * self._K_p * e_star - - Parameters - ---------- - stress : array-like float - The stress - load : array-like float - The load - - Returns - ------- - strain : array-like float - The resulting strain - ''' - - if not isinstance(stress, float): - stress = stress.astype(float) - - return self._ramberg_osgood_relation.strain(stress) - - def load(self, stress, rtol=1e-4, tol=1e-4): - '''Apply the notch-approximation law "backwards", i.e., compute the linear-elastic stress (called "load" or "L" in FKM nonlinear) - from the elastic-plastic stress as from the notch approximation. - This backward step is needed for the pfp FKM nonlinear surface layer & roughness. - - This method is the inverse operation of "stress", i.e., L = load(stress(L)) and S = stress(load(stress)). - - Note that this method is only implemented for the scalar case, as the FKM nonlinear surface layer & roughness - also only handles the scalar case with one assessment point at once, not with entire meshes. - - Parameters - ---------- - stress : array-like float - The elastic-plastic stress as computed by the notch approximation - - Returns - ------- - load : array-like float - The resulting load or linear elastic stress. - - ''' - - x0 = stress / (1 - (1 - 1/self._K_p)/1000) - - # suppress the divergence warnings - with warnings.catch_warnings(): - warnings.simplefilter("ignore") - - load = optimize.newton(func=self._load_implicit, x0=x0, \ - args=([stress]), rtol=rtol, tol=tol, maxiter=50) - - return load - - def stress_secondary_branch(self, delta_load, rtol=1e-4, tol=1e-4): - '''Calculate the stress on secondary branches in the stress-strain diagram at a given - elastic-plastic stress (load), from a FE computation. - This is done by solving for the root of f(sigma) in eq. 2.8-43 of FKM nonlinear. - - Parameters - ---------- - delta_load : array-like float - The load increment of the hysteresis - - Returns - ------- - delta_stress : array-like float - The resulting stress increment within the hysteresis - - Todo - ---- - - In the future, we can evaluate the runtime performance and try a Newton method instead - of the currently used secant method to speed up the computation. - - .. code:: - - fprime=self._d_stress_secondary_implicit_numeric - - ''' - - # initial value as given by correction document to FKM nonlinear - x0 = delta_load * (1 - (1 - 1/self._K_p)/1000) - - # suppress the divergence warnings - with warnings.catch_warnings(): - warnings.simplefilter("ignore") - - delta_stress = optimize.newton(func=self._stress_secondary_implicit, x0=x0, \ - args=([delta_load]), rtol=rtol, tol=tol, full_output=True, maxiter=50) - - # only for multiple points at once, if some points diverged - if hasattr(delta_stress[1], "iloc") and sum(delta_stress[1]) < len(delta_stress[1]): - delta_stress = self._stress_secondary_fix_not_converged_values(delta_stress, delta_load, x0, rtol, tol) - - return delta_stress[0] - + def _stress_secondary_fix_not_converged_values(self, delta_stress, delta_load, x0, rtol, tol): '''For the values that did not converge in the previous vectorized call to optimize.newton, call optimize.newton again on the scalar value. This usually finds the correct solution.''' @@ -412,59 +466,4 @@ def _stress_secondary_fix_not_converged_values(self, delta_stress, delta_load, x if result[1].converged: delta_stress[0].iloc[index_diverged] = result[0] return delta_stress - - def strain_secondary_branch(self, delta_stress, delta_load): - '''Calculate the strain on secondary branches in the stress-strain diagram at a given stress and load. - The formula is given by eq. 2.8-43 of FKM nonlinear. - - Parameters - ---------- - delta_sigma : array-like float - The stress increment - delta_load : array-like float - The load increment - - Returns - ------- - strain : array-like float - The resulting strain - ''' - - if not isinstance(delta_stress, float): - delta_stress = delta_stress.astype(float) - - return self._ramberg_osgood_relation.delta_strain(delta_stress) - - def load_secondary_branch(self, delta_stress, rtol=1e-4, tol=1e-4): - '''Apply the notch-approximation law "backwards", i.e., compute the linear-elastic stress (called "load" or "L" in FKM nonlinear) - from the elastic-plastic stress as from the notch approximation. - This backward step is needed for the pfp FKM nonlinear surface layer & roughness. - - This method is the inverse operation of "stress", i.e., L = load(stress(L)) and S = stress(load(stress)). - - Note that this method is only implemented for the scalar case, as the FKM nonlinear surface layer & roughness - also only handles the scalar case with one assessment point at once, not with entire meshes. - - Parameters - ---------- - delta_stress : array-like float - The increment of the elastic-plastic stress as computed by the notch approximation - - Returns - ------- - delta_load : array-like float - The resulting load or linear elastic stress. - - ''' - - x0 = delta_stress / (1 - (1 - 1/self._K_p)/1000) - - # suppress the divergence warnings - with warnings.catch_warnings(): - warnings.simplefilter("ignore") - - delta_load = optimize.newton(func=self._load_secondary_implicit, x0=x0, \ - args=([delta_stress]), rtol=rtol, tol=tol, maxiter=20) - - return delta_load \ No newline at end of file diff --git a/src/pylife/strength/damage_parameter.py b/src/pylife/strength/damage_parameter.py index 9454ef21..22d6d095 100644 --- a/src/pylife/strength/damage_parameter.py +++ b/src/pylife/strength/damage_parameter.py @@ -51,6 +51,10 @@ def __init__(self, collective, assessment_parameters): self._compute_values() + @property + def collective(self): + return self._collective + def _compute_values(self): """Compute the P_RAM damage parameter according to FKM nonlinear. """ @@ -78,11 +82,7 @@ def _compute_values(self): # delete temporary columns self._collective.drop(columns = ["discriminant", "k"], inplace=True) - - @property - def collective(self): - return self._collective - + class P_RAJ: def __init__(self, collective, assessment_parameters, component_woehler_curve_P_RAJ): @@ -160,7 +160,10 @@ def __init__(self, collective, assessment_parameters, component_woehler_curve_P_ self._P_RAJ_D_0.reset_index(drop=True, inplace=True) self._compute_values() - + + @property + def collective(self): + def _compute_values(self): """Compute the P_RAJ damage parameter according to FKM nonlinear. """ @@ -262,8 +265,7 @@ def _calculate_fatigue_limit_variables(self, D_akt): #P_RAJ_D = pd.Series(self._P_RAJ_D_0) return a_0, delta_J_eff_th, P_RAJ_D - - + def _compute_crack_opening_loop(self): """compute crack opening strain with history (chapter 2.8.9.3, chapter 2.9.8.1 is better)""" @@ -553,7 +555,5 @@ def optimize(row): self._assessment_parameters.P_RAJ_klass_max = self._calculate_P_RAJ(delta_stress, delta_strain) - - @property - def collective(self): + return self._collective \ No newline at end of file diff --git a/src/pylife/strength/fkm_load_distribution.py b/src/pylife/strength/fkm_load_distribution.py index 622f4270..f5c3abd5 100644 --- a/src/pylife/strength/fkm_load_distribution.py +++ b/src/pylife/strength/fkm_load_distribution.py @@ -94,66 +94,6 @@ class FKMLoadSequence(PylifeSignal): mesh.fkm_load_sequence.scaled_by_constant(2) ''' - def _validate(self): - if len(self._obj) == 0: - raise AttributeError("Load series is empty.") - - def _validate_parameters(self, input_parameters, required_parameters): - - for required_parameter in required_parameters: - if required_parameter not in input_parameters: - raise ValueError(f"Given parameters have to include \"{required_parameter}\".") - - def _get_beta(self, input_parameters): - """ - Compute a scaling factor for assessing a load sequence for a given failure probability, - for details, refer to the FKM nonlinear document. - - The beta factors are also described in "A. Fischer. Bestimmung modifizierter Teilsicherheitsbeiwerte zur semiprobabilistischen - Bemessung von Stahlbetonkonstruktionen im Bestand. TU Kaiserslautern, 2010" - - Parameters - ---------- - input_parameters : pd.Series - The set of assessment parameters, has to contain the assessment failure probability ``input_parameters.P_A``. - This variable has to be one of {1e-7, 1e-6, 1e-5, 7.2e-5, 1e-3, 2.3e-1, 0.5}. - - Raises - ------ - ValueError - If the given failure probability has an invalid value. - - Returns - ------- - float - The value of the beta parameter. - - """ - - if np.isclose(input_parameters.P_A, 1e-7): - return 5.20 - - elif np.isclose(input_parameters.P_A, 1e-6): - return 4.75 - - elif np.isclose(input_parameters.P_A, 1e-5): - return 4.27 - - elif np.isclose(input_parameters.P_A, 7.2e-5): - return 3.8 - - elif np.isclose(input_parameters.P_A, 1e-3): - return 3.09 - - elif np.isclose(input_parameters.P_A, 2.3e-1): - return 0.739 - - elif np.isclose(input_parameters.P_A, 0.5): - return 0 - - else: - raise ValueError(f"P_A={input_parameters.P_A} has to be one of "+"{1e-7, 1e-6, 1e-5, 7.2e-5, 1e-3, 2.3e-1, 0.5}.") - def scaled_by_constant(self, gamma_L): """ Scales the load sequence by the given constant gamma_L. This method basically computes @@ -243,6 +183,66 @@ def maximum_absolute_load(self, max_load_independently_for_nodes=False): return float(L_max) + def _validate(self): + if len(self._obj) == 0: + raise AttributeError("Load series is empty.") + + def _validate_parameters(self, input_parameters, required_parameters): + + for required_parameter in required_parameters: + if required_parameter not in input_parameters: + raise ValueError(f"Given parameters have to include \"{required_parameter}\".") + + def _get_beta(self, input_parameters): + """ + Compute a scaling factor for assessing a load sequence for a given failure probability, + for details, refer to the FKM nonlinear document. + + The beta factors are also described in "A. Fischer. Bestimmung modifizierter Teilsicherheitsbeiwerte zur semiprobabilistischen + Bemessung von Stahlbetonkonstruktionen im Bestand. TU Kaiserslautern, 2010" + + Parameters + ---------- + input_parameters : pd.Series + The set of assessment parameters, has to contain the assessment failure probability ``input_parameters.P_A``. + This variable has to be one of {1e-7, 1e-6, 1e-5, 7.2e-5, 1e-3, 2.3e-1, 0.5}. + + Raises + ------ + ValueError + If the given failure probability has an invalid value. + + Returns + ------- + float + The value of the beta parameter. + + """ + + if np.isclose(input_parameters.P_A, 1e-7): + return 5.20 + + elif np.isclose(input_parameters.P_A, 1e-6): + return 4.75 + + elif np.isclose(input_parameters.P_A, 1e-5): + return 4.27 + + elif np.isclose(input_parameters.P_A, 7.2e-5): + return 3.8 + + elif np.isclose(input_parameters.P_A, 1e-3): + return 3.09 + + elif np.isclose(input_parameters.P_A, 2.3e-1): + return 0.739 + + elif np.isclose(input_parameters.P_A, 0.5): + return 0 + + else: + raise ValueError(f"P_A={input_parameters.P_A} has to be one of "+"{1e-7, 1e-6, 1e-5, 7.2e-5, 1e-3, 2.3e-1, 0.5}.") + @pd.api.extensions.register_dataframe_accessor("fkm_safety_normal_from_stddev") @pd.api.extensions.register_series_accessor("fkm_safety_normal_from_stddev") diff --git a/src/pylife/strength/fkm_nonlinear/assessment_nonlinear_standard.py b/src/pylife/strength/fkm_nonlinear/assessment_nonlinear_standard.py index e46dc911..cc9fe0ca 100644 --- a/src/pylife/strength/fkm_nonlinear/assessment_nonlinear_standard.py +++ b/src/pylife/strength/fkm_nonlinear/assessment_nonlinear_standard.py @@ -30,6 +30,193 @@ (FKM non-linear guideline 2019) ''' +def perform_fkm_nonlinear_assessment(assessment_parameters, load_sequence, calculate_P_RAM=True, calculate_P_RAJ=True): + r"""Perform the lifetime assessment according to FKM nonlinear, using the damage parameters P_RAM and/or P_RAJ. + The assessment can be done for a load sequence on a single point or for multiple points at once, e.g., a FEM mesh. + If multiple points at once are used, it is assumed that the load sequences at all nodes are scaled versions of each + other. + + For an assessment with multiple points at once, the relative stress gradient G can be either specified to be constant + or it can have a different value at every point. + + The FKM nonlinear guideline defines three possible methods to consider the statistical distribution of the load: + + 1. a normal distribution with given standard deviation, $s_L$ + 2. a logarithmic-normal distribution with given standard deviation $LSD_s$ + 3. an unknown distribution, use the constant factor :math:`\gamma_L=1.1` for $P_L = 2.5\%$ + or :math:`\gamma_L=1` for $P_L = 50\%$ or + + If the ``assessment_parameters`̀ contain a value for ``s_L``, the first approach is used (normal distribution). + Else, if the ``assessment_parameters`̀ contain a value for ``LSD_s``, the second approach is used (log-normal distribution). + Else, if only ``P_L`̀ is given a scaling with the according factor is used. The statistical assessment can be skipped + by settings P_A = 0.5 and P_L = 50. + + Parameters + ---------- + assessment_parameters : pandas Series + The parameters that specify the material and the assessment problem. The following parameters are required: + + * ``MatGroupFKM``: string, one of {``Steel``, ``SteelCast``, ``Al_wrought``}. Specifies the considered material group. + * ``FinishingFKM``: string, one of {``none``}, the type of surface finishing (Surface finishing types are not implemented for FKM nonlinear). + * ``R_m``: float [MPa], ultimate tensile strength (de: Zugfestigkeit). + Note that this value can also be estimated from a pre-product nominal value, as described in the FKM document. + * ``K_RP``: float, [-], surface roughness factor, set to 1 for polished surfaces or determine from the given diagrams included in the FKM document. + * ``R_z``: float [um], average roughness (de: mittlere Rauheit), only required if K_RP is not specified directly + * ``P_A``: float. Specifies the failure probability for the assessment (de: auszulegende Ausfallwahrscheinlichkeit). + Note that any value for P_A in (0,1) is possible, not just the fixed values that are defined in the FKM nonlinear + guideline + Set to 0.5 to disable statistical assessment, e.g., to simulate cyclic experiments. + * ``beta``: float, damage index, specify this as an alternative to ``P_A``. + * ``P_L``: float, [%], one of {̀ `2.5``%, ``50``%}, probability of occurence of the specified load sequence + (de: Auftretenswahrscheinlilchkeit der Lastfolge). Usually set to 50 to disable statistical assessment for the + load. + * ``s_L``: float (optional), [MPa] standard deviation of Gaussian distribution for the statistical distribution of the load + * ``LSD_s``: float (optional), [MPa] standard deviation of the lognormal distribution for the statistical distribution of the load + * `̀ c``, float, [MPa/N] factor from reference load with which FE simulation was obtained to computed equivalent stress + (de: Übertragungsfaktor Vergleichsspannung zu Referenzlast im Nachweispunkt) c = sigma_V / L_REF + * ̀ `A_sigma``: float, [mm^2] highly loaded surface area of the component (de: Hochbeanspruchte Oberfläche des Bauteils) + * ``A_ref``: float, [mm^2] (de: Hochbeanspruchte Oberfläche eines Referenzvolumens), usually set to 500 + * ``G``: float, [mm^-1] relative stress gradient (de: bezogener Spannungsgradient). + This value can be either a constant value or a pandas Series with different values for every node. + If a Series is used, the order of the G values in the Series has to match the order of the assessment points in the load sequence. + The actual values of the index are irrelevant. + + Note that the relative stress gradient can be computed as follows: + + .. code:: + + grad = pyLife_mesh.gradient_3D.gradient_of('mises') + + # compute the absolute stress gradient + grad["abs_grad"] = np.linalg.norm(grad, axis = 1) + pylife_mesh = pylife_mesh.join(grad, sort=False) + + # compute scaled stress gradient (bezogener Spannungsgradient) + pylife_mesh["G"] = pylife_mesh.abs_grad / pylife_mesh.mises + + To add the value of G to the ``assessment_parameters``, do the following: + + .. code:: + + # remove element_id + G = pylife_mesh['G'].droplevel("element_id") + + # remove duplicate node entries + G = G[~G.index.duplicated(keep="first")].sort_index() + + assessment_parameters["G"] = G + + * ``K_p``: float, [-] (de: Traglastformzahl) K_p = F_plastic / F_yield (3.1.1). + Note that Seeger-Beste and P_RAJ only work for K_p > 1. + * ``n_bins``: int, optional (default: 200) number of bins or classes for P_RAJ computation. A larger value gives more accurate results but longer runtimes. + load_sequence : pandas Series + A sequential list of loads that are applied to the component. If the assessment should be done for + a single points, this is simply a pandas Series. For multiple points at once, it should be a pandas + DataFrame with a two-level MultiIndex with fields "load_step" and "node_id". + The load_step describes the point in time of the sequence and must be consecutive starting from 0. + The node_id identifies the assessment point or mesh node in every load step. The data frame contains + only one column with the stress at every node. The relation between the loads at every nodes + has to be constant over the load steps, i.e., the load sequences at the nodes are scaled versions + of each other. + + An example is given below: + + .. code:: + + S_v + load_step node_id + 0 1 -51.135208 + 2 28.023306 + 3 30.012435 + 4 -11.698302 + 5 287.099222 + ... ... ... + 14614 287.099222 + 1 1 -51.135208 + ... ... ... + 7 1 -51.135208 + ... ... ... + 14610 -113.355076 + 14611 -43.790024 + 14612 -99.422582 + 14613 -77.195496 + 14614 -90.303717 + + calculate_P_RAM : bool (optional) + Whether to use the P_RAM damage parameter for the assessment. Default: True. + calculate_P_RAJ : bool (optional) + Whether to use the P_RAJ damage parameter for the assessment. Default: True. + + Returns + ------- + result : pandas Series + The asssessment result containing at least the following items: + + * ``P_RAM_is_life_infinite``: (bool) whether we have infinite life (de: Dauerfestigkeit) + * ``P_RAM_lifetime_n_cycles``: (float) lifetime in number of cycles + * ``P_RAM_lifetime_n_times_load_sequence``: (float) lifetime how often the full load sequence can be applied + * ``P_RAJ_is_life_infinite`` (bool) whether we have infinite life (de: Dauerfestigkeit) + * ``P_RAJ_lifetime_n_cycles``: (float) lifetime in number of cycles + * ``P_RAJ_lifetime_n_times_load_sequence``: (float) lifetime how often the full load sequence can be applied + + The result dict contains even more entries which are for further information and debugging purposes, such as + woehler curve objects and collective tables. + + If P_A is set to 0.5 and P_L is set to 50, i.e., no statistical assessment is specified, and if the load sequence + is scalar (i.e., not for an entire FEM mesh), the result contains the following additional values: + + * ``P_RAM_lifetime_N_1ppm``, ``P_RAM_lifetime_N_10``, ``P_RAM_lifetime_N_50̀ `, ``P_RAM_lifetime_N_90``: (float) + lifetimes in numbers of cycles, + for P_A = 1ppm = 1e-6, 10%, 50%, and 90%, according to the assessment defined in the FKM nonlinear + guideline. Note that the guideline does not yield a log-normal distributed lifetime. + Furthermore, the value of ``P_RAM_lifetime_N_50̀ ` is lower than the calculated lifetime + ``P_RAM_lifetime_n_cycles``, because it contains a safety factor even for P_A = 50%. + + * ``P_RAM_N_max_bearable``: (function) A python function + ``N_max_bearable(P_A, clip_gamma=False)`` + that calculates the maximum number of cycles + the component can withstand with the given failure probability. + The parameter ``clip_gamma`` specifies whether the scaling factor gamma_M + will be at least 1.1 (P_RAM) or 1.2 (P_RAJ), as defined + in eq. (2.5-38) (PRAM) / eq. (2.8-38) (PRAJ). + + Note that it holds ``P_RAM_lifetime_N_10`` = P_RAM_N_max_bearable(0.1), + and analogously for the variables for 1ppm, 50%, and 90%. + + * ``P_RAM_failure_probability``: (function) A python function, + ``failure_probability(N)`` that calculates the failure probability for a + given number of cycles. + """ + + # check that gradient G is in the correct format + _assert_G_is_in_correct_format(assessment_parameters) + _check_K_p_is_in_range(assessment_parameters) + + scaled_load_sequence = _scale_load_sequence_according_to_probability(assessment_parameters, load_sequence) + scaled_load_sequence = _scale_load_sequence_by_c_factor(assessment_parameters, scaled_load_sequence) + + assessment_parameters = _calculate_local_parameters(assessment_parameters) + + assessment_parameters, component_woehler_curve_P_RAM, component_woehler_curve_P_RAJ \ + = _compute_component_woehler_curves(assessment_parameters) + + maximum_absolute_load = _get_maximum_absolute_load(assessment_parameters, scaled_load_sequence) + + result = {} + + # HCM rainflow counting and damage computation for P_RAM + if calculate_P_RAM: + result = _compute_lifetimes_P_RAM(assessment_parameters, result, scaled_load_sequence, component_woehler_curve_P_RAM, maximum_absolute_load) + + # HCM rainflow counting and damage computation for P_RAJ + if calculate_P_RAJ: + result = _compute_lifetimes_P_RAJ(assessment_parameters, result, scaled_load_sequence, component_woehler_curve_P_RAJ, maximum_absolute_load) + + # additional quantities + result["assessment_parameters"] = assessment_parameters + + return result + def _assert_G_is_in_correct_format(assessment_parameters): """Check that the related stress gradient G is given in the correct format, either as a single float or as a pandas Series with values for each node @@ -54,7 +241,6 @@ def _check_K_p_is_in_range(assessment_parameters): print("Note, K_p is set to 1 which means only P_RAM can be calculated. " f"To use P_RAJ set K_p > 1, e.g. try K_p = 1.001.") - def _scale_load_sequence_according_to_probability(assessment_parameters, load_sequence): """Scales the given load sequence according to one of three methods defined in the FKM nonlinear guideline. @@ -437,190 +623,3 @@ def _compute_lifetimes_P_RAM(assessment_parameters, result, scaled_load_sequence result = _store_additional_objects_in_result_RAM(result, recorder, damage_calculator, component_woehler_curve_P_RAM, detector, detector_1st) return result - -def perform_fkm_nonlinear_assessment(assessment_parameters, load_sequence, calculate_P_RAM=True, calculate_P_RAJ=True): - r"""Perform the lifetime assessment according to FKM nonlinear, using the damage parameters P_RAM and/or P_RAJ. - The assessment can be done for a load sequence on a single point or for multiple points at once, e.g., a FEM mesh. - If multiple points at once are used, it is assumed that the load sequences at all nodes are scaled versions of each - other. - - For an assessment with multiple points at once, the relative stress gradient G can be either specified to be constant - or it can have a different value at every point. - - The FKM nonlinear guideline defines three possible methods to consider the statistical distribution of the load: - - 1. a normal distribution with given standard deviation, $s_L$ - 2. a logarithmic-normal distribution with given standard deviation $LSD_s$ - 3. an unknown distribution, use the constant factor :math:`\gamma_L=1.1` for $P_L = 2.5\%$ - or :math:`\gamma_L=1` for $P_L = 50\%$ or - - If the ``assessment_parameters`̀ contain a value for ``s_L``, the first approach is used (normal distribution). - Else, if the ``assessment_parameters`̀ contain a value for ``LSD_s``, the second approach is used (log-normal distribution). - Else, if only ``P_L`̀ is given a scaling with the according factor is used. The statistical assessment can be skipped - by settings P_A = 0.5 and P_L = 50. - - Parameters - ---------- - assessment_parameters : pandas Series - The parameters that specify the material and the assessment problem. The following parameters are required: - - * ``MatGroupFKM``: string, one of {``Steel``, ``SteelCast``, ``Al_wrought``}. Specifies the considered material group. - * ``FinishingFKM``: string, one of {``none``}, the type of surface finishing (Surface finishing types are not implemented for FKM nonlinear). - * ``R_m``: float [MPa], ultimate tensile strength (de: Zugfestigkeit). - Note that this value can also be estimated from a pre-product nominal value, as described in the FKM document. - * ``K_RP``: float, [-], surface roughness factor, set to 1 for polished surfaces or determine from the given diagrams included in the FKM document. - * ``R_z``: float [um], average roughness (de: mittlere Rauheit), only required if K_RP is not specified directly - * ``P_A``: float. Specifies the failure probability for the assessment (de: auszulegende Ausfallwahrscheinlichkeit). - Note that any value for P_A in (0,1) is possible, not just the fixed values that are defined in the FKM nonlinear - guideline - Set to 0.5 to disable statistical assessment, e.g., to simulate cyclic experiments. - * ``beta``: float, damage index, specify this as an alternative to ``P_A``. - * ``P_L``: float, [%], one of {̀ `2.5``%, ``50``%}, probability of occurence of the specified load sequence - (de: Auftretenswahrscheinlilchkeit der Lastfolge). Usually set to 50 to disable statistical assessment for the - load. - * ``s_L``: float (optional), [MPa] standard deviation of Gaussian distribution for the statistical distribution of the load - * ``LSD_s``: float (optional), [MPa] standard deviation of the lognormal distribution for the statistical distribution of the load - * `̀ c``, float, [MPa/N] factor from reference load with which FE simulation was obtained to computed equivalent stress - (de: Übertragungsfaktor Vergleichsspannung zu Referenzlast im Nachweispunkt) c = sigma_V / L_REF - * ̀ `A_sigma``: float, [mm^2] highly loaded surface area of the component (de: Hochbeanspruchte Oberfläche des Bauteils) - * ``A_ref``: float, [mm^2] (de: Hochbeanspruchte Oberfläche eines Referenzvolumens), usually set to 500 - * ``G``: float, [mm^-1] relative stress gradient (de: bezogener Spannungsgradient). - This value can be either a constant value or a pandas Series with different values for every node. - If a Series is used, the order of the G values in the Series has to match the order of the assessment points in the load sequence. - The actual values of the index are irrelevant. - - Note that the relative stress gradient can be computed as follows: - - .. code:: - - grad = pyLife_mesh.gradient_3D.gradient_of('mises') - - # compute the absolute stress gradient - grad["abs_grad"] = np.linalg.norm(grad, axis = 1) - pylife_mesh = pylife_mesh.join(grad, sort=False) - - # compute scaled stress gradient (bezogener Spannungsgradient) - pylife_mesh["G"] = pylife_mesh.abs_grad / pylife_mesh.mises - - To add the value of G to the ``assessment_parameters``, do the following: - - .. code:: - - # remove element_id - G = pylife_mesh['G'].droplevel("element_id") - - # remove duplicate node entries - G = G[~G.index.duplicated(keep="first")].sort_index() - - assessment_parameters["G"] = G - - * ``K_p``: float, [-] (de: Traglastformzahl) K_p = F_plastic / F_yield (3.1.1). - Note that Seeger-Beste and P_RAJ only work for K_p > 1. - * ``n_bins``: int, optional (default: 200) number of bins or classes for P_RAJ computation. A larger value gives more accurate results but longer runtimes. - load_sequence : pandas Series - A sequential list of loads that are applied to the component. If the assessment should be done for - a single points, this is simply a pandas Series. For multiple points at once, it should be a pandas - DataFrame with a two-level MultiIndex with fields "load_step" and "node_id". - The load_step describes the point in time of the sequence and must be consecutive starting from 0. - The node_id identifies the assessment point or mesh node in every load step. The data frame contains - only one column with the stress at every node. The relation between the loads at every nodes - has to be constant over the load steps, i.e., the load sequences at the nodes are scaled versions - of each other. - - An example is given below: - - .. code:: - - S_v - load_step node_id - 0 1 -51.135208 - 2 28.023306 - 3 30.012435 - 4 -11.698302 - 5 287.099222 - ... ... ... - 14614 287.099222 - 1 1 -51.135208 - ... ... ... - 7 1 -51.135208 - ... ... ... - 14610 -113.355076 - 14611 -43.790024 - 14612 -99.422582 - 14613 -77.195496 - 14614 -90.303717 - - calculate_P_RAM : bool (optional) - Whether to use the P_RAM damage parameter for the assessment. Default: True. - calculate_P_RAJ : bool (optional) - Whether to use the P_RAJ damage parameter for the assessment. Default: True. - - Returns - ------- - result : pandas Series - The asssessment result containing at least the following items: - - * ``P_RAM_is_life_infinite``: (bool) whether we have infinite life (de: Dauerfestigkeit) - * ``P_RAM_lifetime_n_cycles``: (float) lifetime in number of cycles - * ``P_RAM_lifetime_n_times_load_sequence``: (float) lifetime how often the full load sequence can be applied - * ``P_RAJ_is_life_infinite`` (bool) whether we have infinite life (de: Dauerfestigkeit) - * ``P_RAJ_lifetime_n_cycles``: (float) lifetime in number of cycles - * ``P_RAJ_lifetime_n_times_load_sequence``: (float) lifetime how often the full load sequence can be applied - - The result dict contains even more entries which are for further information and debugging purposes, such as - woehler curve objects and collective tables. - - If P_A is set to 0.5 and P_L is set to 50, i.e., no statistical assessment is specified, and if the load sequence - is scalar (i.e., not for an entire FEM mesh), the result contains the following additional values: - - * ``P_RAM_lifetime_N_1ppm``, ``P_RAM_lifetime_N_10``, ``P_RAM_lifetime_N_50̀ `, ``P_RAM_lifetime_N_90``: (float) - lifetimes in numbers of cycles, - for P_A = 1ppm = 1e-6, 10%, 50%, and 90%, according to the assessment defined in the FKM nonlinear - guideline. Note that the guideline does not yield a log-normal distributed lifetime. - Furthermore, the value of ``P_RAM_lifetime_N_50̀ ` is lower than the calculated lifetime - ``P_RAM_lifetime_n_cycles``, because it contains a safety factor even for P_A = 50%. - - * ``P_RAM_N_max_bearable``: (function) A python function - ``N_max_bearable(P_A, clip_gamma=False)`` - that calculates the maximum number of cycles - the component can withstand with the given failure probability. - The parameter ``clip_gamma`` specifies whether the scaling factor gamma_M - will be at least 1.1 (P_RAM) or 1.2 (P_RAJ), as defined - in eq. (2.5-38) (PRAM) / eq. (2.8-38) (PRAJ). - - Note that it holds ``P_RAM_lifetime_N_10`` = P_RAM_N_max_bearable(0.1), - and analogously for the variables for 1ppm, 50%, and 90%. - - * ``P_RAM_failure_probability``: (function) A python function, - ``failure_probability(N)`` that calculates the failure probability for a - given number of cycles. - """ - - # check that gradient G is in the correct format - _assert_G_is_in_correct_format(assessment_parameters) - _check_K_p_is_in_range(assessment_parameters) - - scaled_load_sequence = _scale_load_sequence_according_to_probability(assessment_parameters, load_sequence) - scaled_load_sequence = _scale_load_sequence_by_c_factor(assessment_parameters, scaled_load_sequence) - - assessment_parameters = _calculate_local_parameters(assessment_parameters) - - assessment_parameters, component_woehler_curve_P_RAM, component_woehler_curve_P_RAJ \ - = _compute_component_woehler_curves(assessment_parameters) - - maximum_absolute_load = _get_maximum_absolute_load(assessment_parameters, scaled_load_sequence) - - result = {} - - # HCM rainflow counting and damage computation for P_RAM - if calculate_P_RAM: - result = _compute_lifetimes_P_RAM(assessment_parameters, result, scaled_load_sequence, component_woehler_curve_P_RAM, maximum_absolute_load) - - # HCM rainflow counting and damage computation for P_RAJ - if calculate_P_RAJ: - result = _compute_lifetimes_P_RAJ(assessment_parameters, result, scaled_load_sequence, component_woehler_curve_P_RAJ, maximum_absolute_load) - - # additional quantities - result["assessment_parameters"] = assessment_parameters - - return result diff --git a/src/pylife/strength/fkm_nonlinear/constants.py b/src/pylife/strength/fkm_nonlinear/constants.py index fc47d97a..59568317 100644 --- a/src/pylife/strength/fkm_nonlinear/constants.py +++ b/src/pylife/strength/fkm_nonlinear/constants.py @@ -113,14 +113,6 @@ } }) -def _set_f_25percent(assessment_parameters, resulting_constants): - """Choose the appropriate value for the safety factor f_25%""" - - resulting_constants["f_25percent_material_woehler_RAM"] \ - = resulting_constants["f_25percent_material_woehler_FKM_nonlinear_RAM"] - resulting_constants["f_25percent_material_woehler_RAJ"] \ - = resulting_constants["f_25percent_material_woehler_FKM_nonlinear_RAJ"] - def get_constants_for_material_group(assessment_parameters): """ Retrieve the constants for one of the three material groups that are defined in FKM nonlinear. @@ -153,3 +145,12 @@ def get_constants_for_material_group(assessment_parameters): _set_f_25percent(assessment_parameters, resulting_constants) return resulting_constants + +def _set_f_25percent(assessment_parameters, resulting_constants): + """Choose the appropriate value for the safety factor f_25%""" + + resulting_constants["f_25percent_material_woehler_RAM"] \ + = resulting_constants["f_25percent_material_woehler_FKM_nonlinear_RAM"] + resulting_constants["f_25percent_material_woehler_RAJ"] \ + = resulting_constants["f_25percent_material_woehler_FKM_nonlinear_RAJ"] + \ No newline at end of file diff --git a/src/pylife/strength/fkm_nonlinear/damage_calculator.py b/src/pylife/strength/fkm_nonlinear/damage_calculator.py index c827f691..3d183a36 100644 --- a/src/pylife/strength/fkm_nonlinear/damage_calculator.py +++ b/src/pylife/strength/fkm_nonlinear/damage_calculator.py @@ -87,72 +87,6 @@ def __init__(self, collective, component_woehler_curve_P_RAM): # compute number of cycles until damage sum is 1 self._n_cycles_until_damage = self._collective["cumulative_damage"].groupby("assessment_point_index").apply(lambda array: np.searchsorted(array, 1)) - def _initialize_collective_index(self): - """Assert that the variable `self._collective` contains the proper index - and columns. If the collective does not contain any MultiIndex (or any index at all), - create the appropriate MultiIndex""" - - # if assessment is done for multiple points at once, work with a multi-indexed data frame - if not isinstance(self._collective.index, pd.MultiIndex): - n_hystereses = len(self._collective) - self._collective.index = pd.MultiIndex.from_product([range(n_hystereses), [0]], names=["hysteresis_index", "assessment_point_index"]) - - # assert that the index contains the two columns "hysteresis_index" and "assessment_point_index" - assert self._collective.index.names == ["hysteresis_index", "assessment_point_index"] - - assert "P_RAM" in self._collective - assert "is_closed_hysteresis" in self._collective - assert "run_index" in self._collective - - # store some statistics about the DataFrame - self._n_hystereses = self._collective.groupby("assessment_point_index")["S_min"].count().values[0] - self._n_hystereses_run_2 = self._collective[self._collective["run_index"]==2].groupby("assessment_point_index")["S_min"].count().values[0] - - def _initialize_P_RAM_Z_index(self): - """Properly initialize P_RAM_Z if it was computed individually for every node because of a stress gradient field G. - In such a case, add the proper multi-index with levels "hysteresis_index", "assessment_point_index" such that - the Series is compatible with self._collective. - """ - # if P_RAM_Z is a series without multi-index - if isinstance(self._P_RAM_Z, pd.Series): - if not isinstance(self._P_RAM_Z.index, pd.MultiIndex): - - n_hystereses = len(self._collective.index.get_level_values("hysteresis_index").unique()) - self._P_RAM_Z = pd.Series( - data = (np.ones([n_hystereses,1]) * np.array([self._P_RAM_Z])).flatten(), - index = self._collective.index) - - def _fill_with_default_for_missing_assessment_points(self, df, default_value): - """Add rows to a series df that are not yet there, such that the result has - a row for every assessment point. Example: - - * input ``df``: - - .. code:: - - assessment_point_index - 0 0.312183 - 2 0.312183 - Name: stddev_log_N, dtype: float64 - - * default value: 5 - * result: - - .. code:: - - assessment_point_index - 0 0.312183 - 1 5.000000 - 2 0.312183 - Name: stddev_log_N, dtype: float64 - """ - assessment_point_index = self._collective.index.get_level_values("assessment_point_index").unique() - series_with_all_rows = pd.Series(np.nan, index=assessment_point_index, name="a") - - result = pd.concat([df, series_with_all_rows],axis=1)[df.name] - result = result.fillna(default_value) - return result - @property def collective(self): return self._collective @@ -299,6 +233,72 @@ def failure_probability(N): return N_max_bearable, failure_probability + def _initialize_collective_index(self): + """Assert that the variable `self._collective` contains the proper index + and columns. If the collective does not contain any MultiIndex (or any index at all), + create the appropriate MultiIndex""" + + # if assessment is done for multiple points at once, work with a multi-indexed data frame + if not isinstance(self._collective.index, pd.MultiIndex): + n_hystereses = len(self._collective) + self._collective.index = pd.MultiIndex.from_product([range(n_hystereses), [0]], names=["hysteresis_index", "assessment_point_index"]) + + # assert that the index contains the two columns "hysteresis_index" and "assessment_point_index" + assert self._collective.index.names == ["hysteresis_index", "assessment_point_index"] + + assert "P_RAM" in self._collective + assert "is_closed_hysteresis" in self._collective + assert "run_index" in self._collective + + # store some statistics about the DataFrame + self._n_hystereses = self._collective.groupby("assessment_point_index")["S_min"].count().values[0] + self._n_hystereses_run_2 = self._collective[self._collective["run_index"]==2].groupby("assessment_point_index")["S_min"].count().values[0] + + def _initialize_P_RAM_Z_index(self): + """Properly initialize P_RAM_Z if it was computed individually for every node because of a stress gradient field G. + In such a case, add the proper multi-index with levels "hysteresis_index", "assessment_point_index" such that + the Series is compatible with self._collective. + """ + # if P_RAM_Z is a series without multi-index + if isinstance(self._P_RAM_Z, pd.Series): + if not isinstance(self._P_RAM_Z.index, pd.MultiIndex): + + n_hystereses = len(self._collective.index.get_level_values("hysteresis_index").unique()) + self._P_RAM_Z = pd.Series( + data = (np.ones([n_hystereses,1]) * np.array([self._P_RAM_Z])).flatten(), + index = self._collective.index) + + def _fill_with_default_for_missing_assessment_points(self, df, default_value): + """Add rows to a series df that are not yet there, such that the result has + a row for every assessment point. Example: + + * input ``df``: + + .. code:: + + assessment_point_index + 0 0.312183 + 2 0.312183 + Name: stddev_log_N, dtype: float64 + + * default value: 5 + * result: + + .. code:: + + assessment_point_index + 0 0.312183 + 1 5.000000 + 2 0.312183 + Name: stddev_log_N, dtype: float64 + """ + assessment_point_index = self._collective.index.get_level_values("assessment_point_index").unique() + series_with_all_rows = pd.Series(np.nan, index=assessment_point_index, name="a") + + result = pd.concat([df, series_with_all_rows],axis=1)[df.name] + result = result.fillna(default_value) + return result + class DamageCalculatorPRAJ: """This class performs the lifetime assessment according to the FKM nonlinear assessment. @@ -377,6 +377,112 @@ def __init__(self, collective, assessment_parameters, component_woehler_curve_P_ self._x_bar = (2 + self._xbar_minus_2) + @property + def collective(self): + return self._collective + + @property + def P_RAJ_max(self): + """The maximum P_RAJ damage parameter value of the second run of the HCM algorithm. + If this value is lower than the fatigue strength limit, the component has infinite life. + The method ``compute_damage`` needs to be called beforehand.""" + + # get maximum damage parameter + if isinstance(self._collective.index, pd.MultiIndex): + P_RAJ_max = self._collective.loc[self._collective["run_index"]==2, "P_RAJ"].groupby("assessment_point_index").max() + else: + P_RAJ_max = self._collective.loc[self._collective["run_index"]==2, "P_RAJ"].max() + + return P_RAJ_max.squeeze() + + @property + def is_life_infinite(self): + """Whether the component has infinite life. + The method ``compute_damage`` needs to be called beforehand.""" + + result = self.P_RAJ_max <= self._component_woehler_curve_P_RAJ.fatigue_strength_limit + return result.squeeze() + + @property + def lifetime_n_times_load_sequence(self): + """The number of times the whole load sequence can be traversed until failure. + The method ``compute_damage`` needs to be called beforehand.""" + + # If damage sum of D=1 is reached before end of second run of HCM algorithm, set lifetime_n_times_load_sequence to 0. + # Else the value is x + 1 + result = np.where(self._n_cycles_until_damage < self._n_hystereses, + 0, + self._x_bar) + return result.squeeze() + + @property + def lifetime_n_cycles(self): + """The number of load cycles (as defined in the load collective) until failure. + The method ``compute_damage`` needs to be called beforehand.""" + + # if damage sum of D=1 is reached before end of second run of HCM algorithm + lifetime = np.where(self._n_cycles_until_damage < self._n_hystereses, + self._n_cycles_until_damage, + self._N_bar) + + return lifetime.squeeze() + + def get_lifetime_functions(self): + """Return two python functions that can be used for detailed probabilistic assessment. + + The first function, ``N_max_bearable(P_A, clip_gamma=False)``, + calculates the maximum number of cycles + that the component can withstand with the given failure probability. + The parameter ``clip_gamma`` specifies whether the scaling factor gamma_M + will be at least 1.1 (P_RAM) or 1.2 (P_RAJ), as defined + in eq. (2.5-38) (PRAM) / eq. (2.8-38) (PRAJ) The second function, ``failure_probability(N)``, calculates the failure probability for a + given number of cycles. + + Parameters + ---------- + assessment_parameters : pd.Series + The assessment parameters, only the material group is required to determine the respective + f_2,5% constant. + + Returns + ------- + N_max_bearable + python function, ``N_max_bearable(P_A)`` + failure_probability + python function, ``failure_probability(N)`` + """ + + constants = pylife.strength.fkm_nonlinear.constants.get_constants_for_material_group(self._assessment_parameters) + f_25 = constants.f_25percent_material_woehler_RAJ + slope_woehler = abs(1/self._component_woehler_curve_P_RAJ.d) + lifetime_n_cycles = self.lifetime_n_cycles + + def N_max_bearable(P_A, clip_gamma=False): + beta = pylife.strength.fkm_nonlinear.parameter_calculations.compute_beta(P_A) + log_gamma_M = (0.8*beta - 2)*0.155 + + # Note that the FKM nonlinear guideline defines a cap at 1.2 for P_RAJ. + if clip_gamma: + log_gamma_M = max(log_gamma_M, np.log10(1.2)) + + reduction_factor_P = np.log10(f_25) - log_gamma_M + reduction_factor_N = reduction_factor_P * slope_woehler + + return lifetime_n_cycles * 10**(reduction_factor_N) + + def failure_probability(N): + + result = scipy.optimize.minimize_scalar( + lambda x: (N_max_bearable(x) - N) ** 2, + bounds=[1e-9, 1-1e-9], method='bounded', options={'xatol': 1e-10}) + + if result.success: + return result.x + else: + return 0 + + return N_max_bearable, failure_probability + def _initialize_collective_index(self): """Assert that the variable `self._collective` contains the proper index and columns. If the collective does not contain any MultiIndex (or any index at all), @@ -575,109 +681,4 @@ def f(j): (f(j+1) - f(j)) / denominator, np.inf), 0) - - @property - def collective(self): - return self._collective - - @property - def P_RAJ_max(self): - """The maximum P_RAJ damage parameter value of the second run of the HCM algorithm. - If this value is lower than the fatigue strength limit, the component has infinite life. - The method ``compute_damage`` needs to be called beforehand.""" - - # get maximum damage parameter - if isinstance(self._collective.index, pd.MultiIndex): - P_RAJ_max = self._collective.loc[self._collective["run_index"]==2, "P_RAJ"].groupby("assessment_point_index").max() - else: - P_RAJ_max = self._collective.loc[self._collective["run_index"]==2, "P_RAJ"].max() - - return P_RAJ_max.squeeze() - - @property - def is_life_infinite(self): - """Whether the component has infinite life. - The method ``compute_damage`` needs to be called beforehand.""" - - result = self.P_RAJ_max <= self._component_woehler_curve_P_RAJ.fatigue_strength_limit - return result.squeeze() - - @property - def lifetime_n_times_load_sequence(self): - """The number of times the whole load sequence can be traversed until failure. - The method ``compute_damage`` needs to be called beforehand.""" - - # If damage sum of D=1 is reached before end of second run of HCM algorithm, set lifetime_n_times_load_sequence to 0. - # Else the value is x + 1 - result = np.where(self._n_cycles_until_damage < self._n_hystereses, - 0, - self._x_bar) - return result.squeeze() - - @property - def lifetime_n_cycles(self): - """The number of load cycles (as defined in the load collective) until failure. - The method ``compute_damage`` needs to be called beforehand.""" - - # if damage sum of D=1 is reached before end of second run of HCM algorithm - lifetime = np.where(self._n_cycles_until_damage < self._n_hystereses, - self._n_cycles_until_damage, - self._N_bar) - - return lifetime.squeeze() - - def get_lifetime_functions(self): - """Return two python functions that can be used for detailed probabilistic assessment. - - The first function, ``N_max_bearable(P_A, clip_gamma=False)``, - calculates the maximum number of cycles - that the component can withstand with the given failure probability. - The parameter ``clip_gamma`` specifies whether the scaling factor gamma_M - will be at least 1.1 (P_RAM) or 1.2 (P_RAJ), as defined - in eq. (2.5-38) (PRAM) / eq. (2.8-38) (PRAJ) The second function, ``failure_probability(N)``, calculates the failure probability for a - given number of cycles. - - Parameters - ---------- - assessment_parameters : pd.Series - The assessment parameters, only the material group is required to determine the respective - f_2,5% constant. - - Returns - ------- - N_max_bearable - python function, ``N_max_bearable(P_A)`` - failure_probability - python function, ``failure_probability(N)`` - """ - - constants = pylife.strength.fkm_nonlinear.constants.get_constants_for_material_group(self._assessment_parameters) - f_25 = constants.f_25percent_material_woehler_RAJ - slope_woehler = abs(1/self._component_woehler_curve_P_RAJ.d) - lifetime_n_cycles = self.lifetime_n_cycles - - def N_max_bearable(P_A, clip_gamma=False): - beta = pylife.strength.fkm_nonlinear.parameter_calculations.compute_beta(P_A) - log_gamma_M = (0.8*beta - 2)*0.155 - - # Note that the FKM nonlinear guideline defines a cap at 1.2 for P_RAJ. - if clip_gamma: - log_gamma_M = max(log_gamma_M, np.log10(1.2)) - - reduction_factor_P = np.log10(f_25) - log_gamma_M - reduction_factor_N = reduction_factor_P * slope_woehler - - return lifetime_n_cycles * 10**(reduction_factor_N) - - def failure_probability(N): - - result = scipy.optimize.minimize_scalar( - lambda x: (N_max_bearable(x) - N) ** 2, - bounds=[1e-9, 1-1e-9], method='bounded', options={'xatol': 1e-10}) - - if result.success: - return result.x - else: - return 0 - - return N_max_bearable, failure_probability + \ No newline at end of file diff --git a/src/pylife/stress/rainflow/fkm_nonlinear.py b/src/pylife/stress/rainflow/fkm_nonlinear.py index 34577de8..ed57c0e0 100644 --- a/src/pylife/stress/rainflow/fkm_nonlinear.py +++ b/src/pylife/stress/rainflow/fkm_nonlinear.py @@ -106,108 +106,6 @@ def __init__(self, recorder, notch_approximation_law): self._strain_values = [] self._n_strain_values_first_run = 0 - - def _proceed_on_primary_branch(self, current_point): - """Follow the primary branch (de: Erstbelastungskurve) of a notch approximation material curve. - - Parameters - ---------- - previous_point : _HCM_Point - The starting point in the stress-strain diagram where to begin to follow the primary branch. - current_point : _HCM_Point - The end point until where to follow the primary branch. This variable only needs to have the load value. - - Returns - ------- - current_point : _HCM_Point - The initially given current point, but with updated values of stress and strain. - - """ - sigma = self._notch_approximation_law.stress(current_point.load) - epsilon = self._notch_approximation_law.strain(sigma, current_point.load) - - current_point._stress = sigma - current_point._strain = epsilon - - - # log point for later plotting - self._hcm_point_history.append(("primary", current_point, self._hysteresis_index)) - - return current_point - - def _proceed_on_secondary_branch(self, previous_point, current_point): - """Follow the secondary branch of a notch approximation material curve. - - Parameters - ---------- - previous_point : _HCM_Point - The starting point in the stress-strain diagram where to begin to follow the primary branch. - current_point : _HCM_Point - The end point until where to follow the primary branch. This variable only needs to have the load value. - - Returns - ------- - current_point : _HCM_Point - The initially given current point, but with updated values of stress and strain. - - """ - delta_L = current_point.load - previous_point.load # as described in FKM nonlinear - - delta_sigma = self._notch_approximation_law.stress_secondary_branch(delta_L) - delta_epsilon = self._notch_approximation_law.strain_secondary_branch(delta_sigma, delta_L) - - current_point._stress = previous_point._stress + delta_sigma - current_point._strain = previous_point._strain + delta_epsilon - - # log point for later plotting - self._hcm_point_history.append(("secondary", current_point, self._hysteresis_index)) - - return current_point - - def _adjust_samples_and_flush_for_hcm_first_run(self, samples): - - # convert from Series to np.array - if isinstance(samples, pd.Series): - samples = samples.to_numpy() - - # prepend zero - if isinstance(samples, np.ndarray): - samples = np.concatenate([[0], samples]) - - elif isinstance(samples, pd.DataFrame): - - # get the index with all node_id`s - node_id_index = samples.groupby("node_id").first().index - - # create a new sample with 0 load for all nodes - multi_index = pd.MultiIndex.from_product([[0], node_id_index], names=["load_step","node_id"]) - first_sample = pd.DataFrame(index=multi_index, data=[0]*len(multi_index), columns=samples.columns) - - # increase the load_step index value by one for all samples - samples_without_index = samples.reset_index() - samples_without_index.load_step += 1 - samples = samples_without_index.set_index(["load_step", "node_id"]) - - # prepend the new zero load sample - samples = pd.concat([first_sample, samples], axis=0) - - # determine first, second and last samples - scalar_samples = samples - - if isinstance(samples, pd.DataFrame): - # convert to list - scalar_samples = [df.iloc[0].values[0] for _,df in samples.groupby("load_step")] - - - scalar_samples_twice = np.concatenate([scalar_samples, scalar_samples]) - turn_indices,_ = pylife.stress.rainflow.general.find_turns(scalar_samples_twice) - - flush = True - if len(scalar_samples)-1 not in turn_indices: - flush = False - - return samples, flush - def process_hcm_first(self, samples): """Perform the HCM algorithm for the first time. This processes the given samples accordingly, only considering @@ -341,6 +239,198 @@ def process(self, samples, flush=False): return self + @property + def strain_values(self): + """ + Get the strain values of the turning points in the stress-strain diagram. + They are needed in the FKM nonlinear roughness & surface layer algorithm, which adds residual stresses in another pass of the HCM algorithm. + + Returns + ------- + list of float + The strain values of the turning points that are visited during the HCM algorithm. + """ + + return np.array(self._strain_values) + + @property + def strain_values_first_run(self): + """ + Get the strain values of the turning points in the stress-strain diagram, for the first run of the HCM algorithm. + They are needed in the FKM nonlinear roughness & surface layer algorithm, which adds residual stresses in another pass of the HCM algorithm. + + Returns + ------- + list of float + The strain values of the turning points that are visited during the first run of the HCM algorithm. + """ + + return np.array(self._strain_values[:self._n_strain_values_first_run]) + + @property + def strain_values_second_run(self): + """ + Get the strain values of the turning points in the stress-strain diagram, for the second and any further run of the HCM algorithm. + They are needed in the FKM nonlinear roughness & surface layer algorithm, which adds residual stresses in another pass of the HCM algorithm. + + Returns + ------- + list of float + The strain values of the turning points that are visited during the second run of the HCM algorithm. + """ + + return np.array(self._strain_values[self._n_strain_values_first_run:]) + + def interpolated_stress_strain_data(self, n_points_per_branch=100, only_hystereses=False): + """Return points on the traversed hysteresis curve, mainly intended for plotting. + The curve including all hystereses, primary and secondary branches is sampled + at a fixed number of points within each hysteresis branch. + These points can be used for plotting. + + The intended use is to generate plots as follows: + + .. code:: python + + fkm_nonlinear_detector.process_hcm_first(...) + sampling_parameter = 100 # choose larger for smoother plot or smaller for lower runtime + strain_values_primary, stress_values_primary, _, strain_values_secondary, stress_values_secondary, _ \ + = fkm_nonlinear_detector.interpolated_stress_strain_data(sampling_parameter) + + plt.plot(strain_values_primary, stress_values_primary, "g-", lw=3) + plt.plot(strain_values_secondary, stress_values_secondary, "b-.", lw=1) + + Parameters + ---------- + n_points_per_branch : int, optional + How many sampling points to use per hysteresis branch, default 100. + A larger value values means smoother curves but longer runtime. + + only_hystereses : bool, optional + Default ``False``. If only graphs of the closed hystereses should be output. + Note that this does not work for hysteresis that have multiple smaller hysterseses included. + + Returns + ------- + strain_values_primary, stress_values_primary, hysteresis_index_primary, \ + strain_values_secondary, stress_values_secondary, hysteresis_index_secondary + Lists of strains and stresses of the points on the stress-strain curve, + separately for primary and secondary branches. The lists contain nan values whenever the + curve of the *same* branch is discontinuous. This allows to plot the entire curve + with different colors for primary and secondary branches. + + The entries in hysteresis_index_primary and hysteresis_index_secondary are the row + indices into the collective DataFrame returned by the recorder. This allows, e.g., + to separate the output of multiple runs of the HCM algorithm or to plot the traversed + paths on the stress-strain diagram for individual steps of the algorithm. + + """ + + assert n_points_per_branch >= 2 + + plotter = FKMNonlinearHysteresisPlotter(self._hcm_point_history, self._ramberg_osgood_relation) + return plotter.interpolated_stress_strain_data(n_points_per_branch, only_hystereses) + + def _proceed_on_primary_branch(self, current_point): + """Follow the primary branch (de: Erstbelastungskurve) of a notch approximation material curve. + + Parameters + ---------- + previous_point : _HCM_Point + The starting point in the stress-strain diagram where to begin to follow the primary branch. + current_point : _HCM_Point + The end point until where to follow the primary branch. This variable only needs to have the load value. + + Returns + ------- + current_point : _HCM_Point + The initially given current point, but with updated values of stress and strain. + + """ + sigma = self._notch_approximation_law.stress(current_point.load) + epsilon = self._notch_approximation_law.strain(sigma, current_point.load) + + current_point._stress = sigma + current_point._strain = epsilon + + + # log point for later plotting + self._hcm_point_history.append(("primary", current_point, self._hysteresis_index)) + + return current_point + + def _proceed_on_secondary_branch(self, previous_point, current_point): + """Follow the secondary branch of a notch approximation material curve. + + Parameters + ---------- + previous_point : _HCM_Point + The starting point in the stress-strain diagram where to begin to follow the primary branch. + current_point : _HCM_Point + The end point until where to follow the primary branch. This variable only needs to have the load value. + + Returns + ------- + current_point : _HCM_Point + The initially given current point, but with updated values of stress and strain. + + """ + delta_L = current_point.load - previous_point.load # as described in FKM nonlinear + + delta_sigma = self._notch_approximation_law.stress_secondary_branch(delta_L) + delta_epsilon = self._notch_approximation_law.strain_secondary_branch(delta_sigma, delta_L) + + current_point._stress = previous_point._stress + delta_sigma + current_point._strain = previous_point._strain + delta_epsilon + + # log point for later plotting + self._hcm_point_history.append(("secondary", current_point, self._hysteresis_index)) + + return current_point + + def _adjust_samples_and_flush_for_hcm_first_run(self, samples): + + # convert from Series to np.array + if isinstance(samples, pd.Series): + samples = samples.to_numpy() + + # prepend zero + if isinstance(samples, np.ndarray): + samples = np.concatenate([[0], samples]) + + elif isinstance(samples, pd.DataFrame): + + # get the index with all node_id`s + node_id_index = samples.groupby("node_id").first().index + + # create a new sample with 0 load for all nodes + multi_index = pd.MultiIndex.from_product([[0], node_id_index], names=["load_step","node_id"]) + first_sample = pd.DataFrame(index=multi_index, data=[0]*len(multi_index), columns=samples.columns) + + # increase the load_step index value by one for all samples + samples_without_index = samples.reset_index() + samples_without_index.load_step += 1 + samples = samples_without_index.set_index(["load_step", "node_id"]) + + # prepend the new zero load sample + samples = pd.concat([first_sample, samples], axis=0) + + # determine first, second and last samples + scalar_samples = samples + + if isinstance(samples, pd.DataFrame): + # convert to list + scalar_samples = [df.iloc[0].values[0] for _,df in samples.groupby("load_step")] + + + scalar_samples_twice = np.concatenate([scalar_samples, scalar_samples]) + turn_indices,_ = pylife.stress.rainflow.general.find_turns(scalar_samples_twice) + + flush = True + if len(scalar_samples)-1 not in turn_indices: + flush = False + + return samples, flush + def _perform_hcm_algorithm(self, samples, recording_lists, largest_point, previous_load, iz, ir, load_max_seen, load_turning_points): """Perform the entire HCM algorithm for all load samples, record the found hysteresis parameters in the recording_lists.""" @@ -645,97 +735,6 @@ def _initialize_epsilon_min_for_hcm_run(self, samples, load_turning_points): self._epsilon_min_LF *= pd.Series(1, index=np.arange(n_nodes)) self._epsilon_max_LF *= pd.Series(1, index=np.arange(n_nodes)) - @property - def strain_values(self): - """ - Get the strain values of the turning points in the stress-strain diagram. - They are needed in the FKM nonlinear roughness & surface layer algorithm, which adds residual stresses in another pass of the HCM algorithm. - - Returns - ------- - list of float - The strain values of the turning points that are visited during the HCM algorithm. - """ - - return np.array(self._strain_values) - - @property - def strain_values_first_run(self): - """ - Get the strain values of the turning points in the stress-strain diagram, for the first run of the HCM algorithm. - They are needed in the FKM nonlinear roughness & surface layer algorithm, which adds residual stresses in another pass of the HCM algorithm. - - Returns - ------- - list of float - The strain values of the turning points that are visited during the first run of the HCM algorithm. - """ - - return np.array(self._strain_values[:self._n_strain_values_first_run]) - - @property - def strain_values_second_run(self): - """ - Get the strain values of the turning points in the stress-strain diagram, for the second and any further run of the HCM algorithm. - They are needed in the FKM nonlinear roughness & surface layer algorithm, which adds residual stresses in another pass of the HCM algorithm. - - Returns - ------- - list of float - The strain values of the turning points that are visited during the second run of the HCM algorithm. - """ - - return np.array(self._strain_values[self._n_strain_values_first_run:]) - - def interpolated_stress_strain_data(self, n_points_per_branch=100, only_hystereses=False): - """Return points on the traversed hysteresis curve, mainly intended for plotting. - The curve including all hystereses, primary and secondary branches is sampled - at a fixed number of points within each hysteresis branch. - These points can be used for plotting. - - The intended use is to generate plots as follows: - - .. code:: python - - fkm_nonlinear_detector.process_hcm_first(...) - sampling_parameter = 100 # choose larger for smoother plot or smaller for lower runtime - strain_values_primary, stress_values_primary, _, strain_values_secondary, stress_values_secondary, _ \ - = fkm_nonlinear_detector.interpolated_stress_strain_data(sampling_parameter) - - plt.plot(strain_values_primary, stress_values_primary, "g-", lw=3) - plt.plot(strain_values_secondary, stress_values_secondary, "b-.", lw=1) - - Parameters - ---------- - n_points_per_branch : int, optional - How many sampling points to use per hysteresis branch, default 100. - A larger value values means smoother curves but longer runtime. - - only_hystereses : bool, optional - Default ``False``. If only graphs of the closed hystereses should be output. - Note that this does not work for hysteresis that have multiple smaller hysterseses included. - - Returns - ------- - strain_values_primary, stress_values_primary, hysteresis_index_primary, \ - strain_values_secondary, stress_values_secondary, hysteresis_index_secondary - Lists of strains and stresses of the points on the stress-strain curve, - separately for primary and secondary branches. The lists contain nan values whenever the - curve of the *same* branch is discontinuous. This allows to plot the entire curve - with different colors for primary and secondary branches. - - The entries in hysteresis_index_primary and hysteresis_index_secondary are the row - indices into the collective DataFrame returned by the recorder. This allows, e.g., - to separate the output of multiple runs of the HCM algorithm or to plot the traversed - paths on the stress-strain diagram for individual steps of the algorithm. - - """ - - assert n_points_per_branch >= 2 - - plotter = FKMNonlinearHysteresisPlotter(self._hcm_point_history, self._ramberg_osgood_relation) - return plotter.interpolated_stress_strain_data(n_points_per_branch, only_hystereses) - class FKMNonlinearHysteresisPlotter: def __init__(self, hcm_point_history, ramberg_osgood_relation): diff --git a/src/pylife/stress/rainflow/recorders.py b/src/pylife/stress/rainflow/recorders.py index 7f57f8d5..36875ae9 100644 --- a/src/pylife/stress/rainflow/recorders.py +++ b/src/pylife/stress/rainflow/recorders.py @@ -220,18 +220,6 @@ def epsilon_m(self): return np.where(self.is_zero_mean_stress_and_strain, \ 0, 0.5 * (np.array(self._epsilon_min) + np.array(self._epsilon_max))) - def _get_for_every_node(self, boolean_array): - - # number of points, i.e., number of values for every load step - m = len(self._S_min[0]) - - # bring the array of boolean values to the right shape - # numeric_array contains only 0s and 1s for False and True - numeric_array = np.array(boolean_array).reshape(-1,1).dot(np.ones((1,m))) - - # transform the array to boolean type - return np.where(numeric_array == 1, True, False) - @property def is_zero_mean_stress_and_strain(self): @@ -374,3 +362,15 @@ def record_values_fkm_nonlinear(self, loads_min, loads_max, S_min, S_max, epsilo self._debug_output += debug_output self._run_index += [run_index] * len(S_min) + + def _get_for_every_node(self, boolean_array): + + # number of points, i.e., number of values for every load step + m = len(self._S_min[0]) + + # bring the array of boolean values to the right shape + # numeric_array contains only 0s and 1s for False and True + numeric_array = np.array(boolean_array).reshape(-1,1).dot(np.ones((1,m))) + + # transform the array to boolean type + return np.where(numeric_array == 1, True, False)