From 9724ce03addbec416dc0cc564392b8912738b032 Mon Sep 17 00:00:00 2001 From: Marcus Filipe Date: Wed, 10 May 2023 15:08:15 -0300 Subject: [PATCH 1/8] adding the solver option --- ross/rotor_assembly.py | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/ross/rotor_assembly.py b/ross/rotor_assembly.py index e19bef816..517b00840 100644 --- a/ross/rotor_assembly.py +++ b/ross/rotor_assembly.py @@ -1846,7 +1846,7 @@ def run_unbalance_response( return forced_response - def time_response(self, speed, F, t, ic=None): + def time_response(self, speed, F, t, solver="space-state", ic=None): """Time response for a rotor. This method returns the time response for a rotor @@ -1880,8 +1880,11 @@ def time_response(self, speed, F, t, ic=None): >>> rotor.time_response(speed, F, t) # doctest: +ELLIPSIS (array([0. , 0.18518519, 0.37037037, ... """ - lti = self._lti(speed) - return signal.lsim(lti, F, t, X0=ic) + if solver in ["space-state"]: + lti = self._lti(speed) + return signal.lsim(lti, F, t, X0=ic) + elif solver in ["modal"]: + return None def plot_rotor(self, nodes=1, check_sld=False, length_units="m", **kwargs): """Plot a rotor object. @@ -2371,7 +2374,7 @@ def run_level1(self, n=5, stiffness_range=None, num=5, **kwargs): return results - def run_time_response(self, speed, F, t): + def run_time_response(self, speed, F, t, solver="space-state"): """Calculate the time response. This function will take a rotor object and calculate its time response @@ -2391,6 +2394,12 @@ def run_time_response(self, speed, F, t): Each column corresponds to a dof and each row to a time. t : array Time array. + solver : str, optional + Select which solver is going to be used. Options are: + + 'space-state' : solution will be calculated using the space-state formulation + + 'modal' : the modal reduction will be applied to find the solution Returns ------- @@ -2420,7 +2429,7 @@ def run_time_response(self, speed, F, t): >>> # plot orbit response - plotting 3D orbits - full rotor model: >>> fig3 = response.plot_3d() """ - t_, yout, xout = self.time_response(speed, F, t) + t_, yout, xout = self.time_response(speed, F, t, solver) results = TimeResponseResults(self, t, yout, xout) From be864aecb38ad94ba4449afaad0d75c447150323 Mon Sep 17 00:00:00 2001 From: Marcus Filipe Date: Wed, 10 May 2023 17:20:09 -0300 Subject: [PATCH 2/8] modal integration running inside ross --- ross/defects/integrate_solver.py | 103 ++++++++++++++++++- ross/rotor_assembly.py | 169 ++++++++++++++++++++++++++++++- 2 files changed, 264 insertions(+), 8 deletions(-) diff --git a/ross/defects/integrate_solver.py b/ross/defects/integrate_solver.py index 408b2c254..09154a89f 100644 --- a/ross/defects/integrate_solver.py +++ b/ross/defects/integrate_solver.py @@ -40,13 +40,17 @@ class Integrator: """ - def __init__(self, x0, y0, x, h, func, print_progress=False): + def __init__( + self, x0, y0, x, h, func, print_progress=False, size=24, tolerance=1e-9 + ): self.x0 = x0 self.y0 = y0 self.x = x self.h = h self.func = func self.print_progress = print_progress + self.size = size + self.tolerance = tolerance def rk4(self): # Runge-Kutta 4th order (RK4) @@ -57,7 +61,7 @@ def rk4(self): # Iterate for number of iterations y = self.y0 - result = np.zeros((24, n + 1)) + result = np.zeros((self.size, n + 1)) result[:, 0] = self.y0 # 4th-order Runge-Kutta @@ -90,7 +94,7 @@ def rk45(self): # Iterate for number of iterations y = self.y0 - result = np.zeros((24, n + 1)) + result = np.zeros((self.size, n + 1)) result[:, 0] = self.y0 for i in range(1, n + 1): @@ -138,3 +142,96 @@ def rk45(self): self.x0 = self.x0 + self.h return result + + def rkf45(self): + # Runge-Kutta-Fehlberg (RKF45) + + # Count number of iterations using step size or + # step height h + n = int((self.x - self.x0) / self.h) + + # Iterate for number of iterations + y = self.y0 + result = np.zeros((self.size, n + 1)) + result[:, 0] = self.y0 + + for i in range(1, n + 1): + if i % 10000 == 0 and self.print_progress: + print(f"Iteration: {i} \n Time: {self.x0}") + + h = self.h + found = False + cont = 1 + + "Apply Runge Kutta Formulas to find next value of y" + while not found: + if cont % 10 == 0 and self.print_progress: + print(f"Searched for time step for {cont} times.") + + k1 = 1 * self.func(self.x0, y, i) + + yp2 = y + k1 * (h / 4) + k2 = 1 * self.func(self.x0 + (h / 4), yp2, i) + + yp3 = y + k1 * (3 * h / 32) + k2 * (9 * h / 32) + k3 = 1 * self.func(self.x0 + (3 * h / 8), yp3, i) + + yp4 = ( + y + + k1 * (1932 * h / 2197) + - k2 * (7200 * h / 2197) + + k3 * (7296 * h / 2197) + ) + k4 = 1 * self.func(self.x0 + (12 * h / 13), yp4, i) + + yp5 = ( + y + + k1 * (439 * h / 216) + - k2 * (8 * h) + + k3 * (3680 * h / 513) + - k4 * (845 * h / 4104) + ) + k5 = 1 * self.func(self.x0 + h, yp5, i) + + yp6 = ( + y + - k1 * (8 * h / 27) + + k2 * (2 * h) + - k3 * (3544 * h / 2565) + + k4 * (1859 * h / 4104) + - k5 * (11 * h / 40) + ) + k6 = 1 * self.func(self.x0 + (1 * h / 2), yp6, i) + + # Update next value of y + y = y + h * ( + 25 * k1 / 216 + 1408 * k3 / 2565 + 2197 * k4 / 4104 - 1 * k5 / 5 + ) + + z = y + h * ( + 16 * k1 / 135 + + 6656 * k3 / 12825 + + 28561 * k4 / 56430 + - 9 * k5 / 50 + + 2 * k6 / 55 + ) + + erro = np.max(np.abs(z - y)) + + if erro == 0: + erro = 1e-12 + + s = 0.840896 * (self.tolerance / erro) ** (1 / 4) + + if s > 0.95: + found = True + + h = 0.9 * s * h + + self.x0 += h + + result[:, i] = np.copy(y) + + # Update next value of x + + return result diff --git a/ross/rotor_assembly.py b/ross/rotor_assembly.py index 517b00840..00fba77cf 100644 --- a/ross/rotor_assembly.py +++ b/ross/rotor_assembly.py @@ -1350,6 +1350,139 @@ def _lti(self, speed, frequency=None): return sys + def _modal_solution(self, speed, F, t, number_modes=12): + self.size = 2 * number_modes + + if isinstance(speed, np.ndarray) or isinstance(speed, list): + initial_speed = speed[0] + final_speed = speed[-1] + else: + initial_speed = speed + final_speed = speed + + initial_time = t[0] + dt = t[1] - t[0] + final_time = t[-1] + + # parameters for the time integration + lambdat = 0.00001 + # Faxial = 0 + # TorqueI = 0 + # TorqueF = 0 + + # pre-processing of auxilary variuables for the time integration + sA = ( + initial_speed * np.exp(-lambdat * final_time) + - final_speed * np.exp(-lambdat * initial_time) + ) / (np.exp(-lambdat * final_time) - np.exp(-lambdat * initial_time)) + sB = (final_speed - initial_speed) / ( + np.exp(-lambdat * final_time) - np.exp(-lambdat * initial_time) + ) + + # This code below here is used for acceleration and torque application to the rotor. As of + # september/2020 it is unused, but might be implemented in future releases. These would be + # run-up and run-down operations and variations of operating conditions. + # + # sAT = ( + # TorqueI * np.exp(-lambdat * final_time) - TorqueF * np.exp(-lambdat * initial_time) + # ) / (np.exp(-lambdat * final_time) - np.exp(-lambdat * initial_time)) + # sBT = (TorqueF - TorqueI) / ( + # np.exp(-lambdat * final_time) - np.exp(-lambdat * initial_time) + # ) + # + # SpeedV = sA + sB * np.exp(-lambdat * t) + # TorqueV = sAT + sBT * np.exp(-lambdat * t) + # AccelV = -lambdat * sB * np.exp(-lambdat * t) + + # Get default matrices + K = self.K(speed) + C = self.C(speed) + G = self.G() + M = self.M() + + if self.number_dof == 4: + Kst = np.zeros((self.ndof, self.ndof)) + elif self.number_dof == 6: + Kst = self.Kst() + + y0 = np.zeros(self.size) + + angular_position = ( + sA * t - (sB / lambdat) * np.exp(-lambdat * t) + (sB / lambdat) + ) + + self.Omega = sA + sB * np.exp(-lambdat * t) + self.AccelV = -lambdat * sB * np.exp(-lambdat * t) + + _, ModMat = la.eigh(K, M, type=1, turbo=False) + ModMat = ModMat[:, :number_modes] + ModMat = ModMat + + # Modal transformations + self.Mmodal = ((ModMat.T).dot(M)).dot(ModMat) + self.Cmodal = ((ModMat.T).dot(C)).dot(ModMat) + self.Gmodal = ((ModMat.T).dot(G)).dot(ModMat) + self.Kmodal = ((ModMat.T).dot(K)).dot(ModMat) + self.Kstmodal = ((ModMat.T).dot(Kst)).dot(ModMat) + + self.Fmodal = (ModMat.T).dot(F.T) + + self.inv_Mmodal = np.linalg.pinv(self.Mmodal) + + from ross.defects import Integrator + + x = Integrator( + x0=initial_time, + y0=y0, + x=final_time, + h=dt, + func=self._equation_of_movement, + size=self.size, + ) + x = x.rk45() + + yout = ModMat.dot(x[: int(self.size / 2), :]) + + xout = np.zeros_like(yout) + + return t, yout.T, xout.T + + def _equation_of_movement(self, T, Y, i): + """Calculates the displacement and velocity using state-space representation in the modal domain. + + Parameters + ---------- + T : float + Iteration time. + Y : array + Array of displacement and velocity, in the modal domain. + i : int + Iteration step. + + Returns + ------- + new_Y : array + Array of the new displacement and velocity, in the modal domain. + """ + + positions = Y[: int(self.size / 2)] + velocity = Y[int(self.size / 2) :] # velocity ign space state + + # proper equation of movement to be integrated in time + new_V_dot = ( + +self.Fmodal[:, i] + - ((self.Cmodal + self.Gmodal * self.Omega[i])).dot(velocity) + - ((self.Kmodal + self.Kstmodal * self.AccelV[i]).dot(positions)) + ).dot(self.inv_Mmodal) + + new_X_dot = velocity + + new_Y = np.zeros(self.size) + new_Y[: int(self.size / 2)] = new_X_dot + new_Y[int(self.size / 2) :] = new_V_dot + + return new_Y + def transfer_matrix(self, speed=None, frequency=None, modes=None): """Calculate the fer matrix for the frequency response function (FRF). @@ -1846,7 +1979,7 @@ def run_unbalance_response( return forced_response - def time_response(self, speed, F, t, solver="space-state", ic=None): + def time_response(self, speed, F, t, ic=None, solver="space-state", **kwargs): """Time response for a rotor. This method returns the time response for a rotor @@ -1858,8 +1991,21 @@ def time_response(self, speed, F, t, solver="space-state", ic=None): Force array (needs to have the same length as time array). t : array Time array. (must have the same length than lti.B matrix) + solver : str, optional + Select which solver is going to be used. Options are: + + 'space-state' : solution will be calculated using the space-state formulation + + 'modal' : the modal reduction will be applied to find the solution. + ic : array, optional The initial conditions on the state vector (zero by default). + kwargs : dict, optional + If 'modal' option is passed in solver, you can set the number of modes to use in the modal reduction. Such as: + + number_modes : int + + If this parameter is not passed the code will use a default value of 12. Returns ------- @@ -1884,7 +2030,13 @@ def time_response(self, speed, F, t, solver="space-state", ic=None): lti = self._lti(speed) return signal.lsim(lti, F, t, X0=ic) elif solver in ["modal"]: - return None + number_modes = ( + kwargs.get("number_modes") + if kwargs.get("number_modes") is not None + else 12 + ) + + return self._modal_solution(speed, F, t, number_modes) def plot_rotor(self, nodes=1, check_sld=False, length_units="m", **kwargs): """Plot a rotor object. @@ -2374,7 +2526,7 @@ def run_level1(self, n=5, stiffness_range=None, num=5, **kwargs): return results - def run_time_response(self, speed, F, t, solver="space-state"): + def run_time_response(self, speed, F, t, solver="space-state", **kwargs): """Calculate the time response. This function will take a rotor object and calculate its time response @@ -2399,7 +2551,14 @@ def run_time_response(self, speed, F, t, solver="space-state"): 'space-state' : solution will be calculated using the space-state formulation - 'modal' : the modal reduction will be applied to find the solution + 'modal' : the modal reduction will be applied to find the solution. + + kwargs : dict, optional + If 'modal' option is passed in solver, you can set the number of modes to use in the modal reduction. Such as: + + number_modes : int + + If this parameter is not passed the code will use a default value of 12. Returns ------- @@ -2429,7 +2588,7 @@ def run_time_response(self, speed, F, t, solver="space-state"): >>> # plot orbit response - plotting 3D orbits - full rotor model: >>> fig3 = response.plot_3d() """ - t_, yout, xout = self.time_response(speed, F, t, solver) + t_, yout, xout = self.time_response(speed, F, t, solver=solver) results = TimeResponseResults(self, t, yout, xout) From 29ac28410e60673327e1e6be5c11229d4cc6dbc7 Mon Sep 17 00:00:00 2001 From: Marcus Filipe Date: Wed, 10 May 2023 17:35:24 -0300 Subject: [PATCH 3/8] few corrections --- ross/rotor_assembly.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ross/rotor_assembly.py b/ross/rotor_assembly.py index 00fba77cf..7da98082d 100644 --- a/ross/rotor_assembly.py +++ b/ross/rotor_assembly.py @@ -1443,7 +1443,7 @@ def _modal_solution(self, speed, F, t, number_modes=12): yout = ModMat.dot(x[: int(self.size / 2), :]) - xout = np.zeros_like(yout) + xout = np.zeros((2 * self.ndof, len(t))) return t, yout.T, xout.T From eb09c3e449ae6fdc96e359c250f577ad055fb929 Mon Sep 17 00:00:00 2001 From: Marcus Filipe Date: Thu, 11 May 2023 16:05:58 -0300 Subject: [PATCH 4/8] added test to the time_response method and docs --- ross/rotor_assembly.py | 84 ++++++++++++++++++++++++++++--- ross/tests/test_rotor_assembly.py | 76 ++++++++++++++++++++++++++++ 2 files changed, 152 insertions(+), 8 deletions(-) diff --git a/ross/rotor_assembly.py b/ross/rotor_assembly.py index 7da98082d..0a9a291ef 100644 --- a/ross/rotor_assembly.py +++ b/ross/rotor_assembly.py @@ -1350,7 +1350,47 @@ def _lti(self, speed, frequency=None): return sys - def _modal_solution(self, speed, F, t, number_modes=12): + @check_units + def _modal_analysis(self, speed, F, t, number_modes=12, method="rk45"): + """Modal Analysis + + This method is used to perform the time integration of the rotor using the modal reduction. It uses an Runge-Kutta method within ROSS, + + Parameters + ---------- + speed : float, pint.Quantity + Rotor speed. Default unit is rad/s. + F : array + Force array (needs to have the same number of rows as time array). + Each column corresponds to a dof and each row to a time. + t : array + Time array. + number_modes : int + Number of modes considered in the modal reduction. Default is 12. + method : str + Which algorithm to use. Options are: + + 'rk4' : Runge-Kutta, 4th order + + 'rk45' : Runge-Kutta Cash-Karp, 5th order + + 'rkf45' : Runge-Kutta-Fehlberg, 4th and 5th order + + Default is 'rk45'. + + Returns + ------- + t : array + Time values for the output. + yout : array + System response. + xout : array + Time evolution of the state vector. + """ + + if number_modes > self.ndof: + number_modes = self.ndof + self.size = 2 * number_modes if isinstance(speed, np.ndarray) or isinstance(speed, list): @@ -1439,7 +1479,14 @@ def _modal_solution(self, speed, F, t, number_modes=12): func=self._equation_of_movement, size=self.size, ) - x = x.rk45() + + if method in ["rk4", "rk45", "rkf45"]: + x = getattr(x, method)() + else: + warnings.warn( + f"The method did not match any of {['rk4', 'rk45', 'rkf45']}. Using 'rk45'." + ) + x = x.rk45() yout = ModMat.dot(x[: int(self.size / 2), :]) @@ -2001,11 +2048,20 @@ def time_response(self, speed, F, t, ic=None, solver="space-state", **kwargs): ic : array, optional The initial conditions on the state vector (zero by default). kwargs : dict, optional - If 'modal' option is passed in solver, you can set the number of modes to use in the modal reduction. Such as: + If 'modal' option is passed in solver, you can pass the following: number_modes : int + Number of modes considered in the modal reduction. Default is 12. + method : str + Which algorithm to use. Options are: + + 'rk4' : Runge-Kutta, 4th order + + 'rk45' : Runge-Kutta Cash-Karp, 5th order + + 'rkf45' : Runge-Kutta-Fehlberg, 4th and 5th order - If this parameter is not passed the code will use a default value of 12. + Default is 'rk45'. Returns ------- @@ -2035,8 +2091,11 @@ def time_response(self, speed, F, t, ic=None, solver="space-state", **kwargs): if kwargs.get("number_modes") is not None else 12 ) + method = ( + kwargs.get("method") if kwargs.get("method") is not None else "rk45" + ) - return self._modal_solution(speed, F, t, number_modes) + return self._modal_analysis(speed, F, t, number_modes, method) def plot_rotor(self, nodes=1, check_sld=False, length_units="m", **kwargs): """Plot a rotor object. @@ -2554,11 +2613,20 @@ def run_time_response(self, speed, F, t, solver="space-state", **kwargs): 'modal' : the modal reduction will be applied to find the solution. kwargs : dict, optional - If 'modal' option is passed in solver, you can set the number of modes to use in the modal reduction. Such as: + If 'modal' option is passed in solver, you can pass the following: number_modes : int + Number of modes considered in the modal reduction. Default is 12. + method : str + Which algorithm to use. Options are: + + 'rk4' : Runge-Kutta, 4th order + + 'rk45' : Runge-Kutta Cash-Karp, 5th order + + 'rkf45' : Runge-Kutta-Fehlberg, 4th and 5th order - If this parameter is not passed the code will use a default value of 12. + Default is 'rk45'. Returns ------- @@ -2588,7 +2656,7 @@ def run_time_response(self, speed, F, t, solver="space-state", **kwargs): >>> # plot orbit response - plotting 3D orbits - full rotor model: >>> fig3 = response.plot_3d() """ - t_, yout, xout = self.time_response(speed, F, t, solver=solver) + t_, yout, xout = self.time_response(speed, F, t, solver=solver, **kwargs) results = TimeResponseResults(self, t, yout, xout) diff --git a/ross/tests/test_rotor_assembly.py b/ross/tests/test_rotor_assembly.py index e0db857d9..6b472f1f8 100644 --- a/ross/tests/test_rotor_assembly.py +++ b/ross/tests/test_rotor_assembly.py @@ -1732,6 +1732,82 @@ def test_distinct_dof_elements_error(): Rotor(shaft_elem, [disk0, disk1], [bearing0, bearing1], n_eigen=36) +def test_time_response(rotor7): + ndof = rotor7.ndof + node = 3 + speed = 210 + dt = 0.001 + + t = np.arange(0, 0.05 + dt, dt) + F = np.zeros((len(t), ndof)) + + force = [10, 20] + for f in force: + F[:, 4 * node] = f * np.cos(2 * t) + F[:, 4 * node + 1] = f * np.sin(2 * t) + + yout_ss = rotor7.time_response(speed, F, t, solver="space-state")[1] + yout_m = rotor7.time_response(speed, F, t, solver="modal")[1] + + expected_yout = np.array( + [ + 3.58975900e-08, + -6.07916027e-07, + -1.20146761e-06, + -1.04630257e-06, + -5.72170524e-07, + 3.64535467e-07, + 1.30121807e-06, + 2.16817746e-06, + 3.04976236e-06, + 3.70343595e-06, + 4.28103662e-06, + 4.85911656e-06, + 5.49173576e-06, + 6.14018091e-06, + 6.90159286e-06, + 7.73906607e-06, + 8.45176482e-06, + 9.24584112e-06, + 9.98807794e-06, + 1.06500080e-05, + 1.14250893e-05, + 1.21521896e-05, + 1.29387408e-05, + 1.37671245e-05, + 1.45416063e-05, + 1.53836597e-05, + 1.61687207e-05, + 1.68934233e-05, + 1.75985086e-05, + 1.82057604e-05, + 1.87370653e-05, + 1.91820482e-05, + 1.95439473e-05, + 1.98101499e-05, + 1.99801718e-05, + 2.00656631e-05, + 2.00381359e-05, + 1.99315557e-05, + 1.97285413e-05, + 1.94379441e-05, + 1.90955854e-05, + 1.86707406e-05, + 1.82067469e-05, + 1.77043780e-05, + 1.71477468e-05, + 1.65736507e-05, + 1.59597371e-05, + 1.53147058e-05, + 1.46492473e-05, + 1.39508025e-05, + ] + ) + + assert_almost_equal(yout_ss[1:, 0], expected_yout, 6) + assert_almost_equal(yout_m[1:, 0], expected_yout, 6) + + @pytest.fixture def rotor_6dof(): i_d = 0 From 463cb33e38ce75d4d6d8b56bd4d628647d66a022 Mon Sep 17 00:00:00 2001 From: Marcus Filipe Date: Fri, 12 May 2023 11:42:51 -0300 Subject: [PATCH 5/8] automatic setup of number of modes --- ross/rotor_assembly.py | 66 ++++++++++++++++++++++++++++++++---------- 1 file changed, 50 insertions(+), 16 deletions(-) diff --git a/ross/rotor_assembly.py b/ross/rotor_assembly.py index 0a9a291ef..7c95ac4e0 100644 --- a/ross/rotor_assembly.py +++ b/ross/rotor_assembly.py @@ -1351,7 +1351,7 @@ def _lti(self, speed, frequency=None): return sys @check_units - def _modal_analysis(self, speed, F, t, number_modes=12, method="rk45"): + def _modal_solution(self, speed, F, t, number_modes="force", method="rk45"): """Modal Analysis This method is used to perform the time integration of the rotor using the modal reduction. It uses an Runge-Kutta method within ROSS, @@ -1365,9 +1365,13 @@ def _modal_analysis(self, speed, F, t, number_modes=12, method="rk45"): Each column corresponds to a dof and each row to a time. t : array Time array. - number_modes : int - Number of modes considered in the modal reduction. Default is 12. - method : str + number_modes : str, int, optional + Number of modes considered in the modal reduction. It is possible to pass the following string: + + 'force' : the code will set the number of modes based on the rotor natural frequencies. + + Additionally, if it a int is passed the code will use the number of modes defined by the user. Default is 'force'. + method : str, optional Which algorithm to use. Options are: 'rk4' : Runge-Kutta, 4th order @@ -1388,11 +1392,6 @@ def _modal_analysis(self, speed, F, t, number_modes=12, method="rk45"): Time evolution of the state vector. """ - if number_modes > self.ndof: - number_modes = self.ndof - - self.size = 2 * number_modes - if isinstance(speed, np.ndarray) or isinstance(speed, list): initial_speed = speed[0] final_speed = speed[-1] @@ -1445,8 +1444,6 @@ def _modal_analysis(self, speed, F, t, number_modes=12, method="rk45"): elif self.number_dof == 6: Kst = self.Kst() - y0 = np.zeros(self.size) - angular_position = ( sA * t - (sB / lambdat) * np.exp(-lambdat * t) + (sB / lambdat) ) @@ -1454,7 +1451,39 @@ def _modal_analysis(self, speed, F, t, number_modes=12, method="rk45"): self.Omega = sA + sB * np.exp(-lambdat * t) self.AccelV = -lambdat * sB * np.exp(-lambdat * t) - _, ModMat = la.eigh(K, M, type=1, turbo=False) + lambd, ModMat = la.eigh(K, M, type=1, turbo=False) + + wn = np.sqrt(lambd) + n_harmonics = 5 + + up_freq = np.where(wn >= (n_harmonics * speed))[0].tolist() + max_freq = up_freq[0] + + if max_freq % 2: + forced_modes = int(max_freq - 1) + else: + forced_modes = int(max_freq) + + if isinstance(number_modes, str): + if number_modes in ["force"]: + number_modes = forced_modes + else: + raise ValueError( + f"The string passed in the number of modes of does not match : 'force'. \nThe value was: {number_modes}" + ) + else: + if number_modes > self.ndof: + number_modes = self.ndof + + if number_modes < forced_modes: + warnings.warn( + f"The number of modes of {number_modes} is lower than the forced one. This could lead to wrong results." + ) + + self.size = 2 * number_modes + + y0 = np.zeros(self.size) + ModMat = ModMat[:, :number_modes] ModMat = ModMat @@ -2050,8 +2079,13 @@ def time_response(self, speed, F, t, ic=None, solver="space-state", **kwargs): kwargs : dict, optional If 'modal' option is passed in solver, you can pass the following: - number_modes : int - Number of modes considered in the modal reduction. Default is 12. + number_modes : str, int, optional + Number of modes considered in the modal reduction. It is possible to pass the following string: + + 'force' : the code will set the number of modes based on the rotor natural frequencies. + + Additionally, if it a int is passed the code will use the number of modes defined by the user. Default is 'force'. + method : str Which algorithm to use. Options are: @@ -2089,13 +2123,13 @@ def time_response(self, speed, F, t, ic=None, solver="space-state", **kwargs): number_modes = ( kwargs.get("number_modes") if kwargs.get("number_modes") is not None - else 12 + else "force" ) method = ( kwargs.get("method") if kwargs.get("method") is not None else "rk45" ) - return self._modal_analysis(speed, F, t, number_modes, method) + return self._modal_solution(speed, F, t, number_modes, method) def plot_rotor(self, nodes=1, check_sld=False, length_units="m", **kwargs): """Plot a rotor object. From 8792c324a0dd5faf25cf54f95a419b899e9ffb25 Mon Sep 17 00:00:00 2001 From: Marcus Filipe Date: Tue, 16 May 2023 16:44:21 -0300 Subject: [PATCH 6/8] matmul operator applied --- ross/rotor_assembly.py | 59 ++++++++++++++++++++++++++++-------------- 1 file changed, 39 insertions(+), 20 deletions(-) diff --git a/ross/rotor_assembly.py b/ross/rotor_assembly.py index 7c95ac4e0..970b8457b 100644 --- a/ross/rotor_assembly.py +++ b/ross/rotor_assembly.py @@ -1488,32 +1488,49 @@ def _modal_solution(self, speed, F, t, number_modes="force", method="rk45"): ModMat = ModMat # Modal transformations - self.Mmodal = ((ModMat.T).dot(M)).dot(ModMat) - self.Cmodal = ((ModMat.T).dot(C)).dot(ModMat) - self.Gmodal = ((ModMat.T).dot(G)).dot(ModMat) - self.Kmodal = ((ModMat.T).dot(K)).dot(ModMat) - self.Kstmodal = ((ModMat.T).dot(Kst)).dot(ModMat) + self.Mmodal = (ModMat.T @ M) @ ModMat + self.Cmodal = (ModMat.T @ C) @ ModMat + self.Gmodal = (ModMat.T @ G) @ ModMat + self.Kmodal = (ModMat.T @ K) @ ModMat + self.Kstmodal = (ModMat.T @ Kst) @ ModMat - self.Fmodal = (ModMat.T).dot(F.T) + self.Fmodal = ModMat.T @ F.T self.inv_Mmodal = np.linalg.pinv(self.Mmodal) - from ross.defects import Integrator - - x = Integrator( - x0=initial_time, - y0=y0, - x=final_time, - h=dt, - func=self._equation_of_movement, - size=self.size, - ) + self.i = 0 if method in ["rk4", "rk45", "rkf45"]: + from ross.defects import Integrator + + x = Integrator( + x0=initial_time, + y0=y0, + x=final_time, + h=dt, + func=self._equation_of_movement, + size=self.size, + ) x = getattr(x, method)() + + elif method in ["odeint"]: + from scipy.integrate import odeint + + x = odeint(self._equation_of_movement, y0, t, args=(self.i,), tfirst=True) + else: warnings.warn( - f"The method did not match any of {['rk4', 'rk45', 'rkf45']}. Using 'rk45'." + f"The method did not match any of {['rk4', 'rk45', 'rkf45', 'odeint']}. Using 'rk45'." + ) + from ross.defects import Integrator + + x = Integrator( + x0=initial_time, + y0=y0, + x=final_time, + h=dt, + func=self._equation_of_movement, + size=self.size, ) x = x.rk45() @@ -1547,9 +1564,9 @@ def _equation_of_movement(self, T, Y, i): # proper equation of movement to be integrated in time new_V_dot = ( +self.Fmodal[:, i] - - ((self.Cmodal + self.Gmodal * self.Omega[i])).dot(velocity) - - ((self.Kmodal + self.Kstmodal * self.AccelV[i]).dot(positions)) - ).dot(self.inv_Mmodal) + - ((self.Cmodal + self.Gmodal * self.Omega[i]) @ velocity) + - ((self.Kmodal + self.Kstmodal * self.AccelV[i]) @ positions) + ) @ self.inv_Mmodal new_X_dot = velocity @@ -1557,6 +1574,8 @@ def _equation_of_movement(self, T, Y, i): new_Y[: int(self.size / 2)] = new_X_dot new_Y[int(self.size / 2) :] = new_V_dot + self.i += 1 + return new_Y def transfer_matrix(self, speed=None, frequency=None, modes=None): From e28e76f0a74115537807668ab27d74f1e0b9da2f Mon Sep 17 00:00:00 2001 From: Marcus Filipe Date: Tue, 16 May 2023 17:24:05 -0300 Subject: [PATCH 7/8] force changed to auto --- ross/rotor_assembly.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/ross/rotor_assembly.py b/ross/rotor_assembly.py index 970b8457b..562b4b063 100644 --- a/ross/rotor_assembly.py +++ b/ross/rotor_assembly.py @@ -1351,7 +1351,7 @@ def _lti(self, speed, frequency=None): return sys @check_units - def _modal_solution(self, speed, F, t, number_modes="force", method="rk45"): + def _modal_solution(self, speed, F, t, number_modes="auto", method="rk45"): """Modal Analysis This method is used to perform the time integration of the rotor using the modal reduction. It uses an Runge-Kutta method within ROSS, @@ -1465,7 +1465,7 @@ def _modal_solution(self, speed, F, t, number_modes="force", method="rk45"): forced_modes = int(max_freq) if isinstance(number_modes, str): - if number_modes in ["force"]: + if number_modes in ["auto"]: number_modes = forced_modes else: raise ValueError( @@ -2142,7 +2142,7 @@ def time_response(self, speed, F, t, ic=None, solver="space-state", **kwargs): number_modes = ( kwargs.get("number_modes") if kwargs.get("number_modes") is not None - else "force" + else "auto" ) method = ( kwargs.get("method") if kwargs.get("method") is not None else "rk45" From 454aefcdf59d9d87dcd588b94e23e7cd7c373606 Mon Sep 17 00:00:00 2001 From: Marcus Filipe Date: Wed, 17 May 2023 08:31:29 -0300 Subject: [PATCH 8/8] last force changed to auto --- ross/rotor_assembly.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ross/rotor_assembly.py b/ross/rotor_assembly.py index 562b4b063..f14f4278b 100644 --- a/ross/rotor_assembly.py +++ b/ross/rotor_assembly.py @@ -1368,7 +1368,7 @@ def _modal_solution(self, speed, F, t, number_modes="auto", method="rk45"): number_modes : str, int, optional Number of modes considered in the modal reduction. It is possible to pass the following string: - 'force' : the code will set the number of modes based on the rotor natural frequencies. + 'auto' : the code will set the number of modes based on the rotor natural frequencies. Additionally, if it a int is passed the code will use the number of modes defined by the user. Default is 'force'. method : str, optional