Coverage for neuber_correction\neuber_correction.py: 98%
168 statements
« prev ^ index » next coverage.py v7.9.2, created at 2025-08-14 08:34 +0200
1# coding: utf-8
2"""
3This module contains the NeuberCorrection class, which is used to correct the
4stress using the Neuber correction.
5"""
6import math
7from dataclasses import dataclass
8from typing import List
10import matplotlib.pyplot as plt
11import numpy as np
14@dataclass(frozen=True)
15class MaterialForNeuberCorrection:
16 """
17 Class to store the material properties for the Neuber correction.
18 """
20 yield_strength: float
21 sigma_u: float
22 elastic_mod: float
23 eps_u: float
24 yield_offset: float = 0.002
26 def __post_init__(self):
27 """Validate material properties after initialization."""
28 if self.yield_strength <= 0:
29 raise ValueError("yield_strength (yield strength) must be positive")
30 if self.sigma_u <= 0:
31 raise ValueError("sigma_u (tensile strength) must be positive")
32 if self.elastic_mod <= 0:
33 raise ValueError("elastic_mod (Young's modulus) must be positive")
34 if self.eps_u <= 0:
35 raise ValueError("eps_u (strain at UTS) must be positive")
36 if self.sigma_u <= self.yield_strength:
37 raise ValueError(
38 "sigma_u (tensile strength) must be greater than "
39 "yield_strength (yield strength)"
40 )
42 def __hash__(self):
43 """Custom hash method for memoization."""
44 return hash(
45 (
46 self.yield_strength,
47 self.sigma_u,
48 self.elastic_mod,
49 self.eps_u,
50 self.yield_offset,
51 )
52 )
55@dataclass(frozen=True)
56class NeuberSolverSettings:
57 """
58 Class to store the settings for the Neuber correction.
59 """
61 tolerance: float = 1e-6
62 max_iterations: int = 10000
63 memoization_precision: float = 1e-6
65 def __post_init__(self):
66 """Validate settings after initialization."""
67 if self.tolerance <= 0:
68 raise ValueError("tolerance must be positive")
69 if self.max_iterations <= 0:
70 raise ValueError("max_iterations must be positive")
72 def __hash__(self):
73 """Custom hash method for memoization."""
74 return hash((self.tolerance, self.max_iterations, self.memoization_precision))
77class NeuberCorrection:
78 """
79 Class to correct the stress using the Neuber correction.
80 It is a iterative process, so we need to correct the stress until the
81 difference is less than the tolerance
83 Args:
84 yield_strength: The yield strength of the material
85 sigma_u: The tensile strength of the material
86 elastic_mod: The Young's modulus of the material
87 eps_u: The strain at the tensile strength
88 yield_offset: The offset of the yield strength, usually 0.002
89 tolerance: The tolerance of the correction
90 max_iterations: The maximum number of iterations
91 """
93 instances = {}
95 @classmethod
96 def clear_all_instances(cls):
97 """Clear all cached instances and their memoization tables."""
98 cls.instances.clear()
100 def __new__(
101 cls,
102 material: MaterialForNeuberCorrection,
103 settings: NeuberSolverSettings = NeuberSolverSettings(),
104 ):
105 # Check if instance already exists
106 instance_hash = (hash(material), hash(settings))
107 if instance_hash in cls.instances:
108 return cls.instances[instance_hash]
110 # Create new instance
111 instance = super().__new__(cls)
112 instance.hash = instance_hash
113 return instance
115 def __init__(
116 self,
117 material: MaterialForNeuberCorrection,
118 settings: NeuberSolverSettings = NeuberSolverSettings(),
119 ):
120 # Skip initialization if already initialized
121 if hasattr(self, "material"):
122 return
124 self.material = material
125 self.settings = settings
126 self.memoization_table = {} # Keep as dict for O(1) insertion
127 self.memoization_keys = [] # Keep sorted keys for binary search
128 self.__class__.instances[self.hash] = self
130 def _find_cached_stress(self, target_stress: float) -> float | None:
131 """
132 Find a cached stress value within memoization_precision using binary search.
133 Returns the cached result if found, None otherwise.
134 """
135 if not self.memoization_keys:
136 return None
138 # Binary search for the closest stress value
139 left, right = 0, len(self.memoization_keys) - 1
141 while left <= right:
142 mid = (left + right) // 2
143 cached_stress = self.memoization_keys[mid]
145 if abs(target_stress - cached_stress) < self.settings.memoization_precision:
146 return self.memoization_table[cached_stress]
147 if cached_stress < target_stress:
148 left = mid + 1
149 else:
150 right = mid - 1
152 return None
154 def _insert_sorted(self, stress: float, result: float):
155 """
156 Insert stress value into sorted keys list while maintaining order.
157 """
158 # Find insertion point using binary search
159 left, right = 0, len(self.memoization_keys)
161 while left < right:
162 mid = (left + right) // 2
163 if self.memoization_keys[mid] < stress:
164 left = mid + 1
165 else:
166 right = mid
168 # Insert at the correct position
169 self.memoization_keys.insert(left, stress)
170 self.memoization_table[stress] = result
172 def _calculate_ramberg_osgood_parameter_n(self) -> float:
173 """Calculate the Ramberg-Osgood parameter n."""
174 elastic_ultimate = self.material.sigma_u / self.material.elastic_mod
175 plastic_ultimate = self.material.eps_u - elastic_ultimate
176 return (math.log(plastic_ultimate) - math.log(0.002)) / math.log(
177 self.material.sigma_u / self.material.yield_strength
178 )
180 def _calculate_neuber_correction(self, stress: float) -> float:
181 """
182 Calculates the Neuber correction with Ramberg-Osgood equation.
183 Uses Newton-Raphson iteration to find the corrected stress.
184 """
186 # Check memoization table with efficient binary search
187 cached_result = self._find_cached_stress(stress)
188 if cached_result is not None:
189 return cached_result
191 # Calculate Ramberg-Osgood parameter n
192 n = self._calculate_ramberg_osgood_parameter_n()
194 stress_corr = stress
195 final_difference = float("inf")
197 for _ in range(self.settings.max_iterations):
198 # Safeguard against zero or negative stress_corr
199 if stress_corr <= 0:
200 # If stress_corr becomes zero or negative, reset to a small positive value
201 stress_corr = max(stress * 0.1, 1.0)
202 continue
204 # Calculate total strain (elastic + plastic) - full Ramberg-Osgood curve
205 elastic_strain = stress_corr / self.material.elastic_mod
206 plastic_strain = (
207 self.material.yield_offset
208 * (stress_corr / self.material.yield_strength) ** n
209 )
210 total_strain = elastic_strain + plastic_strain
212 # Calculate Neuber strain
213 neuber_strain = (stress**2) / (self.material.elastic_mod * stress_corr)
215 # Check convergence
216 difference = total_strain - neuber_strain
217 abs_difference = abs(difference)
219 if abs_difference < self.settings.tolerance:
220 final_difference = abs_difference
221 break
223 # Calculate derivatives for Newton-Raphson
224 d_elastic = 1 / self.material.elastic_mod
225 d_plastic = (
226 self.material.yield_offset
227 * n
228 * (stress_corr / self.material.yield_strength) ** (n - 1)
229 / self.material.yield_strength
230 )
231 d_total = d_elastic + d_plastic
233 # Calculate d_neuber
234 d_neuber = -(stress**2) / (self.material.elastic_mod * stress_corr**2)
236 # Newton-Raphson update
237 derivative = d_total - d_neuber
238 if abs(derivative) > 1e-10:
239 stress_corr = stress_corr - difference / derivative
240 else:
241 # Fallback to simple bisection
242 stress_corr = stress_corr * (0.99 if difference > 0 else 1.01)
244 # Always store result in memoization table
245 if final_difference < self.settings.tolerance:
246 # Converged successfully
247 self._insert_sorted(stress, stress_corr)
248 return stress_corr
250 raise ValueError(f"Neuber correction failed for stress {stress}")
252 def correct_stress_values(self, stress_values: List[float]) -> List[float]:
253 """
254 Calculates the Neuber correction for the given stress.
255 """
256 return [self._calculate_neuber_correction(stress) for stress in stress_values]
258 def plot_neuber_diagram(
259 self,
260 stress_value: float,
261 show_plot: bool = True,
262 plot_file: str | None = None,
263 plot_pretty_name: str | None = None,
264 ):
265 """
266 Plots the Neuber diagram showing the Neuber hyperbolic curve and
267 Ramberg-Osgood stress-strain curve for the given stress value.
269 Args:
270 stress_value: The elastic stress value to analyze
271 show_plot: Whether to display the plot (default: True)
272 plot_file: Path to the file to save the plot (default: None)
273 plot_pretty_name: The pretty name of the plot (default: None)
274 """
276 plot_pretty_name = f": {plot_pretty_name}" if plot_pretty_name else ""
278 # Calculate Ramberg-Osgood parameter n
279 n = self._calculate_ramberg_osgood_parameter_n()
281 # Calculate the corrected stress
282 corrected_stress = self._calculate_neuber_correction(stress_value)
284 # Create stress range for plotting
285 stress_range = np.linspace(
286 1, max(stress_value * 1.2, self.material.sigma_u), 1000
287 )
289 # Calculate Ramberg-Osgood strain for each stress
290 strains_ro = []
291 for stress in stress_range:
292 elastic_strain = stress / self.material.elastic_mod
293 plastic_strain = (
294 self.material.yield_offset
295 * (stress / self.material.yield_strength) ** n
296 )
297 strains_ro.append(elastic_strain + plastic_strain)
299 # Calculate Neuber hyperbolic curve
300 # Neuber's rule: σ * ε = σ_elastic^2 / E
301 neuber_constant = stress_value**2 / self.material.elastic_mod
302 strains_neuber = neuber_constant / stress_range
304 # Create the plot
305 fig, ax = plt.subplots(figsize=(12, 8))
307 # Plot Ramberg-Osgood curve
308 ax.plot(
309 strains_ro,
310 stress_range,
311 "b-",
312 linewidth=2,
313 label="Ramberg-Osgood Stress-Strain Curve",
314 )
316 # Plot Hooke's law line (elastic line)
317 hooke_strains = np.linspace(0, self.material.eps_u, 100)
318 hooke_stresses = hooke_strains * self.material.elastic_mod
319 ax.plot(
320 hooke_strains,
321 hooke_stresses,
322 "g-",
323 linewidth=1.5,
324 alpha=0.7,
325 label="Hooke's Law (Elastic Line)",
326 )
328 # Plot Neuber hyperbolic curve
329 ax.plot(
330 strains_neuber,
331 stress_range,
332 "r--",
333 linewidth=2,
334 label="Neuber Hyperbolic Curve",
335 )
337 # Mark the elastic stress point
338 elastic_strain = stress_value / self.material.elastic_mod
339 ax.plot(
340 elastic_strain,
341 stress_value,
342 "go",
343 markersize=10,
344 label=f"Elastic Point (σ={stress_value:.0f} MPa)",
345 )
347 # Mark the corrected stress point
348 corrected_strain = corrected_stress / self.material.elastic_mod
349 plastic_strain = (
350 self.material.yield_offset
351 * (corrected_stress / self.material.yield_strength) ** n
352 )
353 corrected_strain += plastic_strain
355 # Choose marker color based on whether corrected stress is below yield
356 if corrected_stress < self.material.yield_strength:
357 marker_color = "orange" # Orange for below yield
358 marker_label = (
359 f"Corrected Point (σ={corrected_stress:.0f} MPa) - BELOW YIELD"
360 )
361 else:
362 marker_color = "mo" # Magenta for above yield
363 marker_label = f"Corrected Point (σ={corrected_stress:.0f} MPa)"
365 ax.plot(
366 corrected_strain,
367 corrected_stress,
368 marker_color,
369 markersize=10,
370 label=marker_label,
371 )
373 # Mark yield strength line
374 ax.axhline(
375 y=self.material.yield_strength,
376 color="orange",
377 linestyle=":",
378 label=f"Yield Strength ({self.material.yield_strength:.0f} MPa)",
379 )
381 # Mark tensile strength line
382 ax.axhline(
383 y=self.material.sigma_u,
384 color="purple",
385 linestyle=":",
386 label=f"Tensile Strength ({self.material.sigma_u:.0f} MPa)",
387 )
389 # Add intersection point
390 intersection_strain = neuber_constant / corrected_stress
391 ax.plot(
392 intersection_strain,
393 corrected_stress,
394 "ko",
395 markersize=8,
396 label="Intersection Point",
397 )
399 # Add vertical and horizontal lines at intersection point
400 ax.axvline(
401 x=intersection_strain,
402 color="black",
403 linestyle="-",
404 alpha=0.5,
405 linewidth=1,
406 )
407 ax.axhline(
408 y=corrected_stress,
409 color="black",
410 linestyle="-",
411 alpha=0.5,
412 linewidth=1,
413 )
415 # Add value annotations
416 ax.annotate(
417 f"σ = {corrected_stress:.0f} MPa",
418 xy=(intersection_strain, corrected_stress),
419 xytext=(intersection_strain + 0.005, corrected_stress + 50),
420 fontsize=10,
421 ha="left",
422 va="bottom",
423 bbox={"boxstyle": "round,pad=0.3", "facecolor": "white", "alpha": 0.8},
424 arrowprops={"arrowstyle": "->", "connectionstyle": "arc3,rad=0"},
425 )
427 ax.annotate(
428 f"ε = {intersection_strain:.4f}",
429 xy=(intersection_strain, corrected_stress),
430 xytext=(intersection_strain + 0.005, corrected_stress - 100),
431 fontsize=10,
432 ha="left",
433 va="top",
434 bbox={"boxstyle": "round,pad=0.3", "facecolor": "white", "alpha": 0.8},
435 arrowprops={"arrowstyle": "->", "connectionstyle": "arc3,rad=0"},
436 )
438 # Customize the plot
439 ax.set_xlabel("Strain ε [-]", fontsize=12)
440 ax.set_ylabel("Stress σ [MPa]", fontsize=12)
441 ax.set_title(
442 f"Neuber Diagram{plot_pretty_name}\n"
443 f"Elastic Stress: {stress_value:.0f} MPa → "
444 f"Corrected Stress: {corrected_stress:.0f} MPa",
445 fontsize=14,
446 fontweight="bold",
447 )
448 ax.grid(True, alpha=0.3)
449 ax.legend(loc="upper right", fontsize=10)
451 # Set reasonable axis limits
452 ax.set_xlim(0, self.material.eps_u * 1.1)
453 ax.set_ylim(0, max(stress_range) * 1.05)
455 plt.tight_layout()
457 if plot_file:
458 plt.savefig(plot_file)
460 if show_plot:
461 plt.show()
462 else:
463 plt.close(fig)
465 return fig, ax