Untitled
#3D principal stress determination #Laoise Reynolds #20400934 import math # ---------------------------------------------------------------- # 1. Define the stress invariants based on the given tensor components. # ---------------------------------------------------------------- # Given stress components: s11 = 70.0 s22 = 150.0 s33 = 250.0 s12 = 100.0 s13 = -100.0 s23 = -150.0 # Compute invariants: I1 = s11 + s22 + s33 I2 = s11*s22 + s22*s33 + s33*s11 - (s12**2 + s13**2 + s23**2) I3 = ( s11*(s22*s33 - s23**2) - s12*(s12*s33 - s13*s23) + s13*(s12*s23 - s13*s22) ) print("Stress Invariants:") print("I1 =", I1) print("I2 =", I2) print("I3 =", I3) # ---------------------------------------------------------------- # 2. Define the polynomial and its derivative. # The characteristic polynomial is: # P(lambda) = lambda^3 - I1*lambda^2 + I2*lambda - I3 = 0. # ---------------------------------------------------------------- # Coefficients of the cubic polynomial: [1, -I1, I2, -I3] coeffs = [1.0, -I1, I2, -I3] def poly(x, coeffs): """Evaluate polynomial with given coefficients at x. Coefficients are ordered from highest to constant term.""" result = 0 degree = len(coeffs) - 1 for i, c in enumerate(coeffs): result += c * x**(degree - i) return result def dpoly(x, coeffs): """Evaluate derivative of the polynomial at x.""" result = 0 degree = len(coeffs) - 1 for i, c in enumerate(coeffs[:-1]): # skip constant term power = degree - i result += power * c * x**(power - 1) return result # ---------------------------------------------------------------- # 3. Newton-Raphson function for finding one root. # ---------------------------------------------------------------- def newton_raphson_poly(initial_guess, coeffs, tol=1e-6, max_iter=100): x = initial_guess for i in range(max_iter): f_val = poly(x, coeffs) df_val = dpoly(x, coeffs) if df_val == 0: print("Zero derivative encountered!") break x_new = x - f_val / df_val if abs(x_new - x) < tol: return x_new, i+1 x = x_new return x, max_iter # ---------------------------------------------------------------- # 4. Deflation of the polynomial (synthetic division) # Given a root r, factor out (x - r) from the polynomial. # ---------------------------------------------------------------- def deflate_poly(coeffs, root): n = len(coeffs) new_coeffs = [coeffs[0]] # start with the leading coefficient for i in range(1, n): new_coeffs.append(coeffs[i] + root * new_coeffs[-1]) # The last element is the remainder (should be ~0); remove it. return new_coeffs[:-1] # ---------------------------------------------------------------- # 5. Find one root using Newton-Raphson. # Choose an initial guess; here we use 50. (You may try other values.) # ---------------------------------------------------------------- initial_guess = 50.0 root1, iterations = newton_raphson_poly(initial_guess, coeffs) print("\nNewton-Raphson on cubic:") print("Found root 1: {:.6f} (in {} iterations)".format(root1, iterations)) # ---------------------------------------------------------------- # 6. Deflate the cubic to obtain a quadratic. # ---------------------------------------------------------------- coeffs_quad = deflate_poly(coeffs, root1) a, b, c = coeffs_quad # Quadratic: a*x^2 + b*x + c = 0 # Solve the quadratic equation. discriminant = b**2 - 4*a*c if discriminant >= 0: root2 = (-b + math.sqrt(discriminant)) / (2*a) root3 = (-b - math.sqrt(discriminant)) / (2*a) else: # In our case with a symmetric stress tensor the eigenvalues are real, # but we include the complex solution branch for completeness. root2 = (-b + complex(0, math.sqrt(-discriminant))) / (2*a) root3 = (-b - complex(0, math.sqrt(-discriminant))) / (2*a) # ---------------------------------------------------------------- # 7. Display the principal stresses (eigenvalues). # ---------------------------------------------------------------- # Gather roots and sort them (smallest first): roots = [root1, root2, root3] # For sorting, if any are complex, sort by real part. roots.sort(key=lambda r: r.real if isinstance(r, complex) else r) print("\nPrincipal Stresses (Eigenvalues):") for i, r in enumerate(roots, start=1): if isinstance(r, complex): print("Principal Stress {}: {:.6f} + {:.6f}i".format(i, r.real, r.imag)) else: print("Principal Stress {}: {:.6f}".format(i, r))
Leave a Comment