From 9547f22b2121a9db982fbd6da189638b9c903ed1 Mon Sep 17 00:00:00 2001 From: devfernandoa Date: Mon, 21 Oct 2024 13:27:46 -0300 Subject: [PATCH] Adding structural-analysis.py --- structural-analysis.py | 225 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 225 insertions(+) create mode 100644 structural-analysis.py diff --git a/structural-analysis.py b/structural-analysis.py new file mode 100644 index 00000000..f89e08ab --- /dev/null +++ b/structural-analysis.py @@ -0,0 +1,225 @@ +import numpy as np + +# 1. Receive the N nodes, connectivities + +n = input('Enter the number of nodes: ') +connectivities = [] +nodes = [] +loads = [] + +for i in range(int(n)): + x = input(f'Enter the x coordinate of node {i}: ') + y = input(f'Enter the y coordinate of node {i}: ') + node = [x, y] + nodes.append(node) + connectivity = input('Enter a list of connectivities for this node separated by commas: ').split(', ') + for i in range(len(connectivity)): + if connectivity[i].isdigit(): + connectivity[i] = int(connectivity[i]) + else: + print(f"Error: '{connectivity[i]}' is not a valid integer.") + break + else: + connectivities.append(connectivity) + loadX = int(input(f'Enter the load applied to node {i} in X: ')) + loadY = int(input(f'Enter the load applied to node {i} in Y: ')) + load = [loadX, loadY] + loads.append(load) + print('\n') + +area = input('Enter the cross-sectional area: ') +elasticity = input('Enter the modulus of elasticity: ') +tension = input('Enter the tensile stress: ') +rupture = input('Enter the rupture stress: ') + +print('\n\n\n') + +element_list = [] +element_matrices = {} +for i in range(len(connectivities)): + for j in range(len(connectivities[i])): + nodes_pair = sorted([i, int(connectivities[i][j])]) + if tuple(nodes_pair) not in element_list: + print(f'Element {i} connected to node {int(connectivities[i][j])} with incidence in {i} and {int(connectivities[i][j])}') + element_list.append(tuple(nodes_pair)) + node1 = nodes[i] + node2 = nodes[int(connectivities[i][j])] + x1, y1 = float(node1[0]), float(node1[1]) + x2, y2 = float(node2[0]), float(node2[1]) + L = np.sqrt((x2 - x1)**2 + (y2 - y1)**2) + c = (x2 - x1) / L + s = (y2 - y1) / L + print(f'\nFor element {i} we have: c = {c} and s = {s}') + k = float(area) * float(elasticity) / L + element_matrix = (float(elasticity) * float(area))/ L * np.array([[c**2, c*s, -c**2, -c*s], [c*s, s**2, -c*s, -s**2], [-c**2, -c*s, c**2, c*s], [-c*s, -s**2, c*s, s**2]]) + element_matrices[tuple(nodes_pair)] = element_matrix + +# Calculating the degree of freedom numbering +dof = [] +for i in range(int(n)): + dof.append([2*i, 2*i+1]) + print('Node', i, '->', dof[i]) + +# Calculating degrees of freedom with restriction +# Restriction at node 0 (element 0 and 1), and at the last node (element last node * 2i + 1) +restricted_dof = [] +restricted_dof.append(dof[0][0]) +restricted_dof.append(dof[0][1]) +restricted_dof.append(dof[-1][1]) +print(restricted_dof) + +# Degree of freedom of nodes with load +loads_dict = {} +for i in range(int(n)): + if (loads[i][0] != 0 or loads[i][1] != 0): + loads_dict[dof[i][0]] = loads[i][0] + loads_dict[dof[i][1]] = loads[i][1] + print("For degree of freedom", dof[i][0], "the load is", loads[i][0], 'in the positive x direction') + print("For degree of freedom", dof[i][1], "the load is", loads[i][1], 'in the positive y direction') + +# 2. Pre-Processing +# a. Assembly of element matrices + +# b. Superposition of matrices - Global stiffness of the structure +K = np.zeros((2 * len(nodes), 2 * len(nodes))) +for (i, j), matrix in element_matrices.items(): + # Get the global degree of freedom indices for the nodes of the element + dof_indices = dof[i] + dof[j] + for m in range(4): + for n in range(4): + K[dof_indices[m]][dof_indices[n]] += matrix[m][n] + +print('\nGlobal stiffness matrix:') +print(K) + +# c. Assembly of the global load vector of the structure +F = np.zeros(2 * len(nodes)) +for dof, load in loads_dict.items(): + F[dof] = load + +# d. Application of boundary conditions +for dof in restricted_dof: + K[dof] = np.zeros(2 * len(nodes)) + K[:,dof] = np.zeros(2 * len(nodes)) + K[dof][dof] = 1 + F[dof] = 0 + +# Ensure that the main diagonal of K is not zero +for i in range(len(K)): + if round(K[i][i], 3) == 0: + K[i][i] = 1 + +K_copy = K.copy() + +print(f'\nInverse Matrix:') + +# Print the K matrix +for i in range(len(K)): + for j in range(len(K[i])): + print(f'{K[i][j]:.2f}', end=' ') + print() + +# e. Solving the system of equations +# Ensuring that F has the same size as the solution vector +F = F.reshape((2 * len(nodes), 1)) + +# Now solving the system +displacements = np.linalg.solve(K, F) + +print('\nNodal force matrix (in N):') +print(F) + +print('\nDisplacements at nodes (in meters):') +for i in range(len(nodes)): + print(f'Node {i}: ({displacements[2*i][0]}, {displacements[2*i+1][0]})') + +# f. Calculating the stress in each element +stress_list = [] +print('\nStress in each element (in Pa):') +for i in range(len(element_list)): + node1 = nodes[element_list[i][0]] + node2 = nodes[element_list[i][1]] + + x1, y1 = float(node1[0]), float(node1[1]) + x2, y2 = float(node2[0]), float(node2[1]) + L = np.sqrt((x2 - x1)**2 + (y2 - y1)**2) + c = (x2 - x1) / L + s = (y2 - y1) / L + + displacement1 = displacements[2*element_list[i][0]:2*element_list[i][0]+2] + displacement2 = displacements[2*element_list[i][1]:2*element_list[i][1]+2] + + stress = (float(elasticity) / L) * np.array([-c, -s, c, s]).dot(np.concatenate((displacement1, displacement2))) + stress_list.append(stress) + print(f'Element {i}: {stress[0]}') + +# g. Calculating the reactions at the support nodes +# Calculating reaction forces through stress and area +reactions = {} + +print('\nReaction forces in each element (in N):') +for i in range(len(element_list)): + node1 = nodes[element_list[i][0]] + node2 = nodes[element_list[i][1]] + x1, y1 = float(node1[0]), float(node1[1]) + x2, y2 = float(node2[0]), float(node2[1]) + L = np.sqrt((x2 - x1)**2 + (y2 - y1)**2) + c = (x2 - x1) / L + s = (y2 - y1) / L + displacement1 = displacements[2*element_list[i][0]:2*element_list[i][0]+2] + displacement2 = displacements[2*element_list[i][1]:2*element_list[i][1]+2] + stress = (float(elasticity) / L) * np.array([-c, -s, c, s]).dot(np.concatenate((displacement1, displacement2))) + reactions[i] = stress * float(area) + print(f'Element {i}: {reactions[i]}') + + +# 4 +def gauss_seidel(A, b, max_iter=90, tol=1e-20): + x = np.zeros_like(b) + for it_count in range(max_iter): + x_new = np.zeros_like(x) + for i in range(A.shape[0]): + s1 = np.dot(A[i, :i], x[:i]) + s2 = np.dot(A[i, i + 1:], x[i + 1:]) + x_new[i] = (b[i] - s1 - s2) / A[i, i] + if np.allclose(x, x_new, rtol=tol): + break + x = x_new + return x + +A = [[3, -0.1, -0.2], [0.1, 7, -0.3], [0.3, -0.2, 10]] +B = [7.85, -19.3, 71.4] + +print('\nSolution of the linear equations system with the Gauss-Seidel method:') +print(gauss_seidel(np.array(A), np.array(B))) + +# Comparing the displacement from linalg (already calculated) and gauss_seidel +print('\nComparison between numpy.linalg.solve and Gauss-Seidel:') +print('Displacements (in meters):') +print('Numpy:', displacements) +print('Gauss-Seidel:', gauss_seidel(K_copy, F)) + +# 5. Analysis of stress conditions in each element +print('\nAnalysis of stress conditions in each element:') + +rupture = float(rupture) + +for i in range(len(element_list)): + # a. Failure by stress + if abs(stress_list[i]) > rupture: + print(f'Element {i}: Failure by stress. Stress of {stress_list[i]} is greater than the rupture stress {rupture}.') + + # b. Failure by buckling + # The critical buckling load (Pcr) for a column with fixed ends is given by Pcr = pi^2*EI/(KL)^2 + # where E is the modulus of elasticity, I is the moment of inertia, K is the column coefficient (1 for fixed ends), and L is the length of the column. + # For a bar with a circular cross-section, I = pi*r^4/4, where r is the radius of the cross-section. + # Assuming the cross-sectional area is a circle, we can calculate the radius as r = sqrt(area/pi). + # Assuming the bar is a column with fixed ends, K = 1. + # Therefore, the critical buckling load is Pcr = pi^2*E*(pi*r^4/4)/(L^2) = pi^3*E*r^4/(4*L^2). + # The buckling stress is then Pcr/area = pi^3*E*r^4/(4*L^2*area) = pi^3*E*r^2/(4*L^2) = pi^3*E/(4*L^2*(area/pi)) = pi^2*E/(4*L^2/r^2) = pi^2*E*r^2/(4*L^2). + # The buckling load is the buckling stress times the area, or Pcr = pi^2*E*r^2/(4*L^2)*area = pi^2*E*area/(4*L^2/r^2) = pi^2*E*area*r^2/(4*L^2). + # Therefore, if the stress in the element is negative (compression) and the magnitude of the stress is greater than the buckling stress, the element fails by buckling. + r = np.sqrt(float(area) / np.pi) + Pcr = np.pi**2 * float(elasticity) * float(area) * r**2 / (4 * L**2) + if stress_list[i] < 0 and abs(stress_list[i]) > Pcr: + print(f'Element {i}: Failure by buckling. Magnitude of compressive stress {abs(stress_list[i])} is greater than the critical buckling load {Pcr}.')