Source code for gastop.evaluator

"""evaluator.py
This file is a part of GASTOp
Authors: Amlan Sinha, Cristian Lacey, Daniel Shaw, Paul Kaneelil, Rory Conlin, Susan Redmond
Licensed under GNU GPLv3.
This module implements the Evaluator class.

"""
import numpy as np


[docs]class Evaluator(): """Implements various methods for scoring the truss in different areas. Methods include calculations of mass, factor of safety, deflections, and interference with user specified areas. The class is designed to be instantiated as an Evaluator object which will fully evaluate a Truss object using specified methods and parameters. """
[docs] def __init__(self, struct_solver, mass_solver, interferences_solver, cost_solver, boundary_conditions, properties_dict): """Creates an Evaluator callable object. Once created, the Evaluator can be called on a Truss object to calculate and assign mass, factor of safety, deflections, etc to the truss. Args: struct_solver (str): Name of the method to be used for structural analysis and calculating fos and deflections, as a string. e.g. ``'mat_struct_analysis_DSM'``. mass_solver (str): Name of the method to be used to calculate mass. e.g. ``'mass_basic'``. interferences_solver (str): Name of method to be used to determine interferences. e.g. ``'interferences_ray_tracing'``. boundary_conditions (dict): Dictionary containing: - ``'loads'`` *(ndarray)*: Array of loads applied to the structure. First index corresponds to the node where the load is applied, second index is the force in x,y,z and moment about x,y,z, third index is for multiple loading scenarios. - ``'fixtures'`` *(ndarray)*: Array of flags denoting whether a node is fixed or free. First index corresponds to the node, the second index corresponds to fixing displacements in x,y,z and rotations about x,y,z. The third index corresponds to multiple loading scenarios with different fixtures for each. Values of the array are 0 (free) or 1 (fixed). properties_dict (dict): Dictionary containing beam properties. Entries should be 1D arrays, with length equal to the number of beam options. Each entry in the array is the value of the key property for the specified beam type. Properties include: - ``'OD'``: Outer diameter of the beam, in meters. - ``'ID'``: Inner diameter of the beam, in meters. - ``'elastic_modulus'``: Elastic or Young's modulus of the material, in Pascals. - ``'yield_strength'``: Yield or failure strength of the material, in Pascals. - ``'shear_modulus'``: Shear modulus of the material, in Pascals. - ``'poisson_ratio'``: Poisson ratio of the material, dimensionless. - ``'x_section_area'``: Cross sectional area of the beam, in square meters. - ``'moment_inertia_y'``: Area moment of inertia about beams y axis, in meters^4. - ``'moment_inertia_z'``: Area moment of inertia about beams z axis, in meters^4. - ``'polar_moment_inertia'``: Area moment of inertia about beams polar axis, in meters^4. - ``'dens'``: Density of the material, in kilograms per cubic meter. cost_solver (str): Name of the method to be used to calculate cost. e.g. ``'cost_calc'``. Returns: callable Evaluator object. """ self.struct_solver = getattr(self, struct_solver) self.mass_solver = getattr(self, mass_solver) self.interferences_solver = getattr(self, interferences_solver) self.boundary_conditions = boundary_conditions self.cost_solver = getattr(self, cost_solver) self.properties_dict = properties_dict
[docs] @staticmethod def mat_struct_analysis_DSM(truss, boundary_conditions, properties_dict): """Calculates deflections and stresses using direct stiffness method. Constructs global stiffness matrix from nodes and connections, and computes deflections under each loading scenario. From deflections, calculates internal forces, stresses, and factor of safety in each member under each loading scenario Args: truss (Truss object): Truss to be evaluated. Must have nodes, edges, and properties defined. boundary_conditions (dict): Dictionary containing: - ``'loads'`` *(ndarray)*: Array of loads applied to the structure. First index corresponds to the node where the load is applied, second index is the force in x,y,z and moment about x,y,z, third index is for multiple loading scenarios. - ``'fixtures'`` *(ndarray)*: Array of flags denoting whether a node is fixed or free. First index corresponds to the node, the second index corresponds to fixing displacements in x,y,z and rotations about x,y,z. The third index corresponds to multiple loading scenarios with different fixtures for each. Values of the array are 0 (free) or 1 (fixed). properties_dict (dict): Dictionary containing beam properties. Entries should be 1D arrays, with length equal to the number of beam options. Each entry in the array is the value of the key property for the specified beam type. Properties include: - ``'OD'``: Outer diameter of the beam, in meters. - ``'ID'``: Inner diameter of the beam, in meters. - ``'elastic_modulus'``: Elastic or Young's modulus of the material, in Pascals. - ``'yield_strength'``: Yield or failure strength of the material, in Pascals. - ``'shear_modulus'``: Shear modulus of the material, in Pascals. - ``'poisson_ratio'``: Poisson ratio of the material, dimensionless. - ``'x_section_area'``: Cross sectional area of the beam, in square meters. - ``'moment_inertia_y'``: Area moment of inertia about beams y axis, in meters^4. - ``'moment_inertia_z'``: Area moment of inertia about beams z axis, in meters^4. - ``'polar_moment_inertia'``: Area moment of inertia about beams polar axis, in meters^4. - ``'dens'``: Density of the material, in kilograms per cubic meter. Returns: 2-element tuple containing: - **fos** *(ndarray)*: 2D array of factor of safety values. First index corresponds to members, second index corresponds to different loading scenarios. Factor of safety is defined as the materials yield strength divided by the von Mises stress in the member. If structure is statically indeterminate under a given loading scenario, fos will be zero. Factor of safety in member i under loading j is fos[i, j] - **deflections** *(ndarray)*: 3D array of node deflections. Distances in meters, angles in radians. First index corresponds to node number, second index is deflections in global x,y,z coordinates, and rotations about global x,y,z axes. The third axis corresponds to different loading scenarios. Deflection at node i under loading j is deflections[i, :, j] = [dx, dy, dz, d_theta_x, d_theta_y, d_theta_z] """ nodes, con, matl = truss.cleaned_params() loads = boundary_conditions['loads'].copy() fixtures = boundary_conditions['fixtures'].copy() eps = np.finfo(float).eps # machine precision num_nodes = nodes.shape[0] num_con = con.shape[0] num_loads = loads.shape[2] # get material properties etc E = properties_dict['elastic_modulus'][matl] G = properties_dict['shear_modulus'][matl] YS = properties_dict['yield_strength'][matl] A = properties_dict['x_section_area'][matl] Iz = properties_dict['moment_inertia_z'][matl] Iy = properties_dict['moment_inertia_y'][matl] J = properties_dict['polar_moment_inertia'][matl] OD = properties_dict['outer_diameter'][matl] # initialize empty matrices # member stiffness matrices in local coords Kloc = np.zeros((12, 12, num_con)) # member stiffness matrices in global coords KlocT = np.zeros((12, 12, num_con)) Kglob = np.zeros((6*num_nodes, 6*num_nodes)) # global stiffness matrix Q = np.zeros((num_con, 12)) # end forces on members Ei = np.zeros((num_con, 12)) # local to global matrix indices V = np.zeros((num_nodes, 6, num_loads)) # displacements r = np.zeros((3, 3, num_con)) # member rotation matrices # local to global transformation matrix T = np.zeros((12, 12, num_con)) k1 = np.zeros((3, 3, num_con)) # stiffness matrix components k2 = np.zeros((3, 3, num_con)) # stiffness matrix components k3 = np.zeros((3, 3, num_con)) # stiffness matrix components k4 = np.zeros((3, 3, num_con)) # stiffness matrix components FoS = np.zeros((num_con, num_loads)) # factor of safety # calculate length and direction sines/cosines of members in global coords edge_vec = nodes[con[:, 1], :] - nodes[con[:, 0], :] # length of projection in x-y plane rho = np.sqrt(edge_vec[:, 0]**2+edge_vec[:, 2]**2) L = np.sqrt(rho**2 + edge_vec[:, 1]**2) # total length sing = rho < eps # id vectors along z axis, azimuth is not defined, convention set to 0 edge_vec[:, 0][sing] = 1 # cos(0) = 1 edge_vec[:, 2][sing] = 0 # sin(0) = 0 rho[sing] = 1 ca = edge_vec[:, 0]/rho # cosine of azimuthal angle sa = edge_vec[:, 2]/rho # sine of azimuthal angle rho[sing] = 0 cp = rho/L # cosine of polar elevation angle sp = edge_vec[:, 1]/L # sine of polar elevation angle # populate transformation matrices r[0, 0, :] = cp*ca r[0, 1, :] = sp r[0, 2, :] = sa*cp r[1, 0, :] = -sp*ca r[1, 1, :] = cp r[1, 2, :] = -sp*sa r[2, 0, :] = -sa r[2, 2, :] = ca T[0:3, 0:3, :] = r T[3:6, 3:6, :] = r T[6:9, 6:9, :] = r T[9:12, 9:12, :] = r # stiffness matrix elements in x,y,z,theta co = np.stack((12*np.ones(num_con), 6*L, 4*L**2, 2*L**2), axis=1) x = (E*A)/L # axial stiffness along x y = (((E*Iy)/(L**3))*co.T).T # bending stiffness about y z = (((E*Iz)/(L**3))*co.T).T # bending stiffness about z g = G*J/L # torsional stiffness about x axis # form stacked local stiffness matrices k1[0, 0, :] = x k1[1, 1, :] = z[:, 0] k1[2, 2, :] = y[:, 0] k2[1, 2, :] = z[:, 1] k2[2, 1, :] = -y[:, 1] k3[0, 0, :] = g k3[1, 1, :] = y[:, 2] k3[2, 2, :] = z[:, 2] k4[0, 0, :] = -g k4[1, 1, :] = y[:, 3] k4[2, 2, :] = z[:, 3] k2t = np.transpose(k2, axes=(1, 0, 2)) kr1 = np.concatenate((k1, k2, -k1, k2), axis=1) kr2 = np.concatenate((k2t, k3, -k2t, k4), axis=1) kr3 = np.concatenate((-k1, -k2, k1, -k2), axis=1) kr4 = np.concatenate((k2t, k4, -k2t, k3), axis=1) Kloc = np.concatenate((kr1, kr2, kr3, kr4), axis=0) # construct global stiffness matrix from element matrices for ii in range(num_con): # get member indices to global stiffness matrix e = np.concatenate((np.arange(6*con[ii, 0], 6*con[ii, 0]+6), np.arange(6*con[ii, 1], 6*con[ii, 1]+6)), axis=0) # transform from local to global coords KlocT[:, :, ii] = np.matmul(Kloc[:, :, ii], T[:, :, ii]) # form global stiffness matrix Kglob[np.ix_(e, e)] = (Kglob[np.ix_(e, e)] + np.matmul(T[:, :, ii].T, KlocT[:, :, ii])) Ei[ii, :] = e # save indices for later # fix free node displacements # currently acts weird and gives wrong results if a node is fixed in one direction and free in another # calculate displacements for j in range(num_loads): # set unconnected unloaded nodes to fixed unconnected = np.setdiff1d( range(num_nodes), np.concatenate((con.flatten(), np.nonzero(loads[:, :, j].any(axis=1))[0]))) fixtures[unconnected, :, j] = 1 # get indices of free nodes f = np.nonzero( 1-np.ravel(fixtures[:, :, j])) f = f[0] # get array out of tuple # solve for displacements of free nodes try: np.ravel(V[:, :, j])[f] = np.linalg.solve( Kglob[np.ix_(f, f)], np.ravel(loads[:, :, j])[f]) # if matrix is singular, stop, FoS still all zeros except np.linalg.LinAlgError: return FoS, V # calculate forces and stresses for i in range(num_con): # end forces Q[i, :] = np.matmul( KlocT[:, :, i], np.ravel(V[:, :, j])[Ei[i, :].astype(int)]) # combined moment about y, z M = np.sqrt(Q[i, 4]**2 + Q[i, 5]**2) # axial stress due to bending moment sigmaXbending = M*OD[i]/(2*Iz[i]) # axial stress due to axial forces sigmaXaxial = np.abs(Q[i, 0]/A[i]) # transverse stress due to torsion tauTorsion = Q[i, 3]*OD[i]/(2*J[i]) # transverse stress due to shear tauXY = 2*np.sqrt(Q[i, 1]**2 + Q[i, 2]**2)/A[i] # determine von mises stress sigmaVM = np.amax((np.sqrt((sigmaXbending+sigmaXaxial)**2 + 3*tauTorsion**2), np.sqrt(sigmaXaxial**2 + 3*tauXY**2))) # factor of safety in each beam under each loading condition if sigmaVM > YS[i]/1000: # to avoid div/0 and huge fos FoS[i, j] = YS[i]/sigmaVM else: FoS[i, j] = 1000 return FoS, V
[docs] @staticmethod def mass_basic(truss, properties_dict): """Calculates mass of structure Considers only members, does not account for additional mass due to welds or connection hardware. Args: truss (Truss object): Truss to be evaluated. Must have nodes, edges, and properties defined. properties_dict (dict): Dictionary containing beam properties. Entries should be 1D arrays, with length equal to the number of beam options. Each entry in the array is the value of the key property for the specified beam type. Properties include: - ``'x_section_area'``: Cross sectional area of the beam, in square meters. - ``'dens'``: Density of the material, in kilograms per cubic meter. Returns: mass (float): Mass of the structure in kilograms. """ nodes, con, matl = truss.cleaned_params() # calculate member lengths L = np.sqrt( np.sum((nodes[con[:, 1], :] - nodes[con[:, 0], :])**2, axis=1)) # get material properties # print(properties_dict['x_section_area']) A = properties_dict['x_section_area'][matl] dens = properties_dict['density'][matl] mass = np.sum(A*L*dens) return mass
[docs] @staticmethod def cost_calc(truss, properties_dict): """Calculates cost of structure Considers only members, does not account for additional cost due to welds or connection hardware. Args: truss (Truss object): Truss to be evaluated. Must have nodes, edges, and properties defined. properties_dict (dict): Dictionary containing beam properties. Entries should be 1D arrays, with length equal to the number of beam options. Each entry in the array is the value of the key property for the specified beam type. Properties include: - ``'x_section_area'``: Cross sectional area of the beam, in square meters. - ``'cost'``: Cost of the material, in $ per meter. Returns: cost (float): Cost of the structure in $. """ nodes, con, matl = truss.cleaned_params() # calculate member lengths edge_vec = nodes[con[:, 1], :] - nodes[con[:, 0], :] L = np.sqrt( edge_vec[:, 0]**2 + edge_vec[:, 1]**2 + edge_vec[:, 2]**2) # get material properties cost_per_len = properties_dict['cost'][matl] mass = np.sum(L*cost_per_len) return mass
[docs] @staticmethod def interference_ray_tracing(truss): """Not implemented yet. TODO: method to determine if truss members are crossing into user specified areas. Used when a structure must be designed around something, such as a passenger compartment or other design components. """ return None
[docs] @staticmethod def blank_test(truss, *args, **kwargs): """Blank function used for testing GA when no evaluation needed Args: truss (Truss object): Dummy Truss object, no attributes required. Returns: 2-element tuple of (None, None) """ return None, None
[docs] def __call__(self, truss): """Computes mass, deflections, etc, and stores it in truss object. Used when an Evaluator object has been created with the methods to be used and any necessary parameters. Args: truss (Truss object): truss to be evaluated. Returns: None """ truss.mark_duplicates() truss.fos, truss.deflection = self.struct_solver( truss, self.boundary_conditions, self.properties_dict) truss.mass = self.mass_solver(truss, self.properties_dict) truss.cost = self.cost_solver(truss, self.properties_dict) truss.interference = self.interferences_solver(truss) return truss