diff --git a/docs/references.bib b/docs/references.bib index 4281c45e5..a31e5ba8b 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -720,3 +720,14 @@ @article{genta1988conical journal = {Journal of Sound and Vibration}, doi = {10.1016/0022-460X(88)90342-2} } + +@article{rouch1980dynamic, +author = {Rouch, K. E. and Kao, J. S.}, +title = {Dynamic Reduction in Rotor Dynamics by the Finite Element Method}, +journal = {Journal of Mechanical Design}, +volume = {102}, +number = {2}, +pages = {360-368}, +year = {1980}, +doi = {10.1115/1.3254752} +} diff --git a/ross/disk_element.py b/ross/disk_element.py index ac226ab95..fcfcd9df8 100644 --- a/ross/disk_element.py +++ b/ross/disk_element.py @@ -90,18 +90,21 @@ def __eq__(self, other): True """ false_number = 0 - for i in self.__dict__: - try: - if np.allclose(self.__dict__[i], other.__dict__[i]): - pass - else: - false_number += 1 - - except TypeError: - if self.__dict__[i] == other.__dict__[i]: - pass - else: - false_number += 1 + if "DiskElement" in other.__class__.__name__: + for i in self.__dict__: + try: + if np.allclose(self.__dict__[i], other.__dict__[i]): + pass + else: + false_number += 1 + + except TypeError: + if self.__dict__[i] == other.__dict__[i]: + pass + else: + false_number += 1 + else: + false_number += 1 if false_number == 0: return True diff --git a/ross/rotor_assembly.py b/ross/rotor_assembly.py index 41f517521..0fef81ebb 100644 --- a/ross/rotor_assembly.py +++ b/ross/rotor_assembly.py @@ -583,7 +583,7 @@ def __eq__(self, other): return False @check_units - def run_modal(self, speed, num_modes=12, sparse=True): + def run_modal(self, speed, num_modes=12, sparse=True, synchronous=False): """Run modal analysis. Method to calculate eigenvalues and eigvectors for a given rotor system. @@ -615,6 +615,9 @@ def run_modal(self, speed, num_modes=12, sparse=True): If False, scipy.linalg.eig() is used to calculate all the eigenvalues and eigenvectors. Default is True. + synchronous : bool, optional + If True a synchronous analysis is carried out according to :cite:`rouch1980dynamic`. + Default is False. Returns ------- @@ -638,7 +641,9 @@ def run_modal(self, speed, num_modes=12, sparse=True): >>> mode2 = 1 # Second mode >>> fig = modal.plot_mode_2d(mode2) """ - evalues, evectors = self._eigen(speed, num_modes=num_modes, sparse=sparse) + evalues, evectors = self._eigen( + speed, num_modes=num_modes, sparse=sparse, synchronous=synchronous + ) wn_len = num_modes // 2 wn = (np.absolute(evalues))[:wn_len] wd = (np.imag(evalues))[:wn_len] @@ -887,9 +892,15 @@ def convergence(self, n_eigval=0, err_max=1e-02): return results - def M(self, frequency=None): + def M(self, frequency=None, synchronous=False): """Mass matrix for an instance of a rotor. + Parameters + ---------- + synchronous : bool, optional + If True a synchronous analysis is carried out. + Default is False. + Returns ------- M0 : np.ndarray @@ -913,10 +924,44 @@ def M(self, frequency=None): for elm in self.elements: dofs = list(elm.dof_global_index.values()) - try: - M0[np.ix_(dofs, dofs)] += elm.M(frequency) - except TypeError: - M0[np.ix_(dofs, dofs)] += elm.M() + if synchronous: + if elm in self.shaft_elements: + try: + M = elm.M(frequency) + except TypeError: + M = elm.M() + G = elm.G() + for i in range(8): + if i == 0 or i == 3 or i == 4 or i == 7: + M[i, 0] = M[i, 0] - G[i, 1] + M[i, 3] = M[i, 3] + G[i, 2] + M[i, 4] = M[i, 4] - G[i, 5] + M[i, 7] = M[i, 7] + G[i, 6] + else: + M[i, 1] = M[i, 1] + G[i, 0] + M[i, 2] = M[i, 2] - G[i, 3] + M[i, 5] = M[i, 5] + G[i, 4] + M[i, 6] = M[i, 6] - G[i, 7] + M0[np.ix_(dofs, dofs)] += M + elif elm in self.disk_elements: + try: + M = elm.M(frequency) + except TypeError: + M = elm.M() + G = elm.G() + M[2, 2] = M[2, 2] - G[2, 3] + M[3, 3] = M[3, 3] + G[3, 2] + M0[np.ix_(dofs, dofs)] += M + else: + try: + M0[np.ix_(dofs, dofs)] += elm.M(frequency) + except TypeError: + M0[np.ix_(dofs, dofs)] += elm.M() + else: + try: + M0[np.ix_(dofs, dofs)] += elm.M(frequency) + except TypeError: + M0[np.ix_(dofs, dofs)] += elm.M() return M0 @@ -1053,7 +1098,7 @@ def G(self): return G0 - def A(self, speed=0, frequency=None): + def A(self, speed=0, frequency=None, synchronous=False): """State space matrix for an instance of a rotor. Parameters @@ -1063,6 +1108,9 @@ def A(self, speed=0, frequency=None): Default is 0. frequency : float, optional Excitation frequency. Default is rotor speed. + synchronous : bool, optional + If True a synchronous analysis is carried out. + Default is False. Returns ------- @@ -1089,7 +1137,7 @@ def A(self, speed=0, frequency=None): # fmt: off A = np.vstack( [np.hstack([Z, I]), - np.hstack([la.solve(-self.M(frequency), self.K(frequency)), la.solve(-self.M(frequency), (self.C(frequency) + self.G() * speed))])]) + np.hstack([la.solve(-self.M(frequency, synchronous=synchronous), self.K(frequency)), la.solve(-self.M(frequency,synchronous=synchronous), (self.C(frequency) + self.G() * speed))])]) # fmt: on return A @@ -1240,7 +1288,14 @@ def _index(eigenvalues): @check_units def _eigen( - self, speed, num_modes=12, frequency=None, sorted_=True, A=None, sparse=True + self, + speed, + num_modes=12, + frequency=None, + sorted_=True, + A=None, + sparse=True, + synchronous=False, ): """Calculate eigenvalues and eigenvectors. @@ -1264,6 +1319,9 @@ def _eigen( sparse : bool, optional If sparse, eigenvalues will be calculated with arpack. Default is True. + synchronous : bool, optional + If True a synchronous analysis is carried out. + Default is False. Returns ------- @@ -1280,20 +1338,34 @@ def _eigen( 91.796... """ if A is None: - A = self.A(speed=speed, frequency=frequency) + A = self.A(speed=speed, frequency=frequency, synchronous=synchronous) - if sparse is True: - try: - evalues, evectors = las.eigs( - A, k=num_modes, sigma=0, ncv=2 * num_modes, which="LM", v0=self._v0 - ) - # store v0 as a linear combination of the previously - # calculated eigenvectors to use in the next call to eigs - self._v0 = np.real(sum(evectors.T)) - except las.ArpackError: - evalues, evectors = la.eig(A) - else: + if synchronous: evalues, evectors = la.eig(A) + idx = np.where(np.imag(evalues) != 0)[0] + evalues = evalues[idx] + evectors = evectors[:, idx] + idx = np.where(np.abs(np.real(evalues) / np.imag(evalues)) < 1000)[0] + evalues = evalues[idx] + evectors = evectors[:, idx] + else: + if sparse is True: + try: + evalues, evectors = las.eigs( + A, + k=num_modes, + sigma=0, + ncv=2 * num_modes, + which="LM", + v0=self._v0, + ) + # store v0 as a linear combination of the previously + # calculated eigenvectors to use in the next call to eigs + self._v0 = np.real(sum(evectors.T)) + except las.ArpackError: + evalues, evectors = la.eig(A) + else: + evalues, evectors = la.eig(A) if sorted_ is False: return evalues, evectors @@ -2176,9 +2248,8 @@ def run_ucs( Default is 16. In this case 4 modes are plotted, since for each pair of eigenvalues calculated we have one wn, and we show only the forward mode in the plots. - synchronous : bool - If True a synchronous analysis is carried out and the frequency of - the first forward model will be equal to the speed. + synchronous : bool, optional + If True a synchronous analysis is carried out. Default is False. Returns @@ -2214,32 +2285,10 @@ def run_ucs( for i, k in enumerate(stiffness_log): bearings = [b.__class__(b.n, kxx=k, cxx=0) for b in bearings_elements] rotor = self.__class__(self.shaft_elements, self.disk_elements, bearings) - speed = 0 - if synchronous: - - def wn_diff(x): - """Function to evaluate difference between speed and - natural frequency for the first mode.""" - modal = rotor.run_modal(speed=x, num_modes=num_modes) - # get first forward mode - if modal.whirl_direction()[0] == "Forward": - wn0 = modal.wn[0] - else: - wn0 = modal.wn[1] - - return wn0 - x - - speed = newton(wn_diff, 0) - modal = rotor.run_modal(speed=speed, num_modes=num_modes) - - # if sync, select only forward modes - if synchronous: - rotor_wn[:, i] = modal.wn[modal.whirl_direction() == "Forward"] - # if not sync, with speed=0 whirl direction can be confusing, with - # two close modes being forward or backward, so we select one mode in - # each 2 modes. - else: - rotor_wn[:, i] = modal.wn[::2] + modal = rotor.run_modal( + speed=0, num_modes=num_modes, synchronous=synchronous + ) + rotor_wn[:, i] = modal.wn[::2] bearing0 = bearings_elements[0] diff --git a/ross/tests/test_rotor_assembly.py b/ross/tests/test_rotor_assembly.py index 73ae9f7b6..c9f8380e2 100644 --- a/ross/tests/test_rotor_assembly.py +++ b/ross/tests/test_rotor_assembly.py @@ -1833,70 +1833,579 @@ def rotor8(): def test_ucs_calc_rotor2(rotor2): + exp_rotor_wn = np.array( + [ + [ + 215.37072557, + 283.80994851, + 366.86912206, + 460.53914551, + 555.6339197, + 640.32650515, + 706.0780713, + 751.34310872, + 779.83176945, + 796.72990074, + 806.39721162, + 811.81289519, + 814.81102149, + 816.45987555, + 817.36339085, + 817.85749941, + 818.1274194, + 818.27478247, + 818.35520922, + 818.39909624, + ], + [ + 598.02474114, + 799.46636352, + 1057.70436175, + 1373.74248509, + 1729.64320906, + 2079.77186424, + 2367.99373458, + 2567.69744007, + 2690.52418862, + 2761.4174174, + 2801.13738888, + 2823.0913233, + 2835.14811303, + 2841.74862093, + 2845.35623074, + 2847.32634838, + 2848.40174108, + 2848.98860223, + 2849.30882097, + 2849.4835338, + ], + [ + 3956.224973, + 4062.07552707, + 4249.58400228, + 4573.6000451, + 5112.58991935, + 5963.9057796, + 7227.66815371, + 8977.98586645, + 11209.15347504, + 13742.89203466, + 16175.83274873, + 18056.34394912, + 19282.21914334, + 20006.24626009, + 20414.15343628, + 20639.66981718, + 20763.42435719, + 20831.12561174, + 20868.11080392, + 20888.30264105, + ], + [ + 4965.28982295, + 5025.38045612, + 5136.67510967, + 5343.17672239, + 5723.17881707, + 6398.10091437, + 7508.44360294, + 9150.62680727, + 11311.73143116, + 13795.63516793, + 16186.47408516, + 18081.26060018, + 19332.61770536, + 20072.40777584, + 20489.35934417, + 20719.90875241, + 20846.4293277, + 20915.646125, + 20953.45866009, + 20974.10006581, + ], + ] + ) ucs_results = rotor2.run_ucs() - expected_intersection_points_y = [ - 215.37072557303264, - 215.37072557303264, - 598.024741157381, - 598.024741157381, - 3956.224979518562, - 3956.224979518562, - 4965.289823255782, - 4965.289823255782, - ] - assert_allclose( - ucs_results.intersection_points["y"], expected_intersection_points_y + assert_allclose(ucs_results.wn, exp_rotor_wn, rtol=1e-6) + + exp_rotor_wn = np.array( + [ + [ + 2.15371329e02, + 2.83812357e02, + 3.66877867e02, + 4.60566615e02, + 5.55704705e02, + 6.40471324e02, + 7.06315117e02, + 7.51667171e02, + 7.80222437e02, + 7.97164913e02, + 8.06859321e02, + 8.12290756e02, + 8.15297782e02, + 8.16951585e02, + 8.17857830e02, + 8.18353436e02, + 8.18624175e02, + 8.18771986e02, + 8.18852657e02, + 8.18896678e02, + ], + [ + 1.88802037e03, + 2.28476088e03, + 2.78000154e03, + 3.41117879e03, + 4.23129187e03, + 5.31089177e03, + 6.73193403e03, + 8.56265868e03, + 1.08004453e04, + 1.32860389e04, + 1.56703469e04, + 1.75831951e04, + 1.88890430e04, + 1.96891387e04, + 2.01514788e04, + 2.04110573e04, + 2.05547675e04, + 2.06337750e04, + 2.06770552e04, + 2.07007192e04, + ], + [ + 4.00701928e03, + 4.11399080e03, + 4.30342165e03, + 4.63058499e03, + 5.17436636e03, + 6.03219769e03, + 7.30332537e03, + 9.05967180e03, + 1.12948065e04, + 1.38430517e04, + 1.63370844e04, + 1.83684523e04, + 1.97674949e04, + 2.06278515e04, + 2.11255558e04, + 2.14050391e04, + 2.15597608e04, + 2.16448157e04, + 2.16914058e04, + 2.17168786e04, + ], + [ + 3.75678825e04, + 3.75949522e04, + 3.76446065e04, + 3.77357411e04, + 3.79031768e04, + 3.82113096e04, + 3.87797759e04, + 3.98314192e04, + 4.17770276e04, + 4.53362190e04, + 5.16143414e04, + 6.19874527e04, + 7.79184129e04, + 1.01059654e05, + 1.33577562e05, + 1.78470328e05, + 2.39880642e05, + 3.23483351e05, + 4.37008883e05, + 5.90956828e05, + ], + ] ) - fig = ucs_results.plot() - assert_allclose(fig.data[4]["y"], expected_intersection_points_y) + ucs_results = rotor2.run_ucs(synchronous=True) + assert_allclose(ucs_results.wn, exp_rotor_wn) def test_ucs_calc(rotor8): - exp_stiffness_range = np.array([1000000.0, 1832980.710832, 3359818.286284]) - exp_rotor_wn = np.array([86.658114, 95.660573, 101.868254]) - exp_intersection_points_x = np.array( - [10058123.652648, 10058123.652648, 10363082.398797] + exp_rotor_wn = np.array( + [ + [ + 86.65811435, + 95.66057326, + 101.86825429, + 105.77854216, + 108.09888143, + 109.42658356, + 110.17043717, + 110.58225348, + 110.80874181, + 110.93285112, + 111.00072364, + 111.03780095, + 111.05804338, + 111.06909117, + 111.07511968, + 111.07840898, + 111.0802036, + 111.08118271, + 111.08171688, + 111.0820083, + ], + [ + 274.31285391, + 325.11814102, + 366.19875009, + 394.79384259, + 412.67726755, + 423.17290554, + 429.12464766, + 432.439225, + 434.26762238, + 435.27109643, + 435.82032762, + 436.12049426, + 436.28441024, + 436.37388292, + 436.42270952, + 436.44935147, + 436.46388748, + 436.47181809, + 436.47614484, + 436.47850536, + ], + [ + 716.78631221, + 838.8349461, + 975.91941882, + 1094.44878855, + 1172.5681763, + 1216.43412434, + 1239.90023448, + 1252.41864548, + 1259.14192099, + 1262.77498349, + 1264.74615415, + 1265.81822821, + 1266.4021086, + 1266.72035103, + 1266.89388148, + 1266.98852606, + 1267.04015246, + 1267.06831513, + 1267.08367898, + 1267.09206063, + ], + [ + 1066.1562004, + 1194.71380733, + 1383.21963246, + 1611.2607587, + 1811.78963328, + 1932.88449588, + 1992.93768733, + 2022.35515478, + 2037.28135654, + 2045.0830944, + 2049.23769603, + 2051.47407412, + 2052.68517943, + 2053.34324267, + 2053.70146233, + 2053.89665559, + 2054.00307642, + 2054.06111426, + 2054.09277044, + 2054.11003892, + ], + ] ) - exp_intersection_points_y = np.array([107.542416, 107.542416, 409.451575]) ucs_results = rotor8.run_ucs() - assert_allclose(ucs_results.stiffness_log[:3], exp_stiffness_range) - assert_allclose(ucs_results.wn[0, :3], exp_rotor_wn) - assert_allclose( - ucs_results.intersection_points["x"][:3], exp_intersection_points_x, rtol=1e-3 - ) - assert_allclose( - ucs_results.intersection_points["y"][:3], exp_intersection_points_y, rtol=1e-3 + assert_allclose(ucs_results.wn, exp_rotor_wn) + + exp_rotor_wn = np.array( + [ + [ + 86.90424993, + 96.07439757, + 102.44368355, + 106.4793955, + 108.88395104, + 110.26337332, + 111.03737524, + 111.46625295, + 111.70223863, + 111.83158678, + 111.9023347, + 111.94098589, + 111.96208851, + 111.97360604, + 111.97989096, + 111.98332019, + 111.98519116, + 111.98621193, + 111.98676882, + 111.98707265, + ], + [ + 288.61681387, + 340.2268304, + 381.1079096, + 409.06350367, + 426.3367896, + 436.40207642, + 442.08713715, + 445.24635051, + 446.98699943, + 447.94170533, + 448.46406291, + 448.7494886, + 448.90533871, + 448.99040391, + 449.03682385, + 449.0621522, + 449.07597137, + 449.08351086, + 449.0876242, + 449.08986829, + ], + [ + 925.78185127, + 1119.90966576, + 1392.50590974, + 1750.68004093, + 2182.0883474, + 2636.00694936, + 3028.14253219, + 3302.1410502, + 3467.15137752, + 3559.91633395, + 3610.85969749, + 3638.66057981, + 3653.81442225, + 3662.07522781, + 3666.57962521, + 3669.03626945, + 3670.37627387, + 3671.10725239, + 3671.50602245, + 3671.72356858, + ], + [ + 1115.3680757, + 1270.72738636, + 1513.13125359, + 1856.25476513, + 2290.15538281, + 2763.19618667, + 3184.21829187, + 3484.87824241, + 3668.02539974, + 3771.46132472, + 3828.3605264, + 3859.43124229, + 3876.37175328, + 3885.60756844, + 3890.64388364, + 3893.39070372, + 3894.8890089, + 3895.7063473, + 3896.15223098, + 3896.39548013, + ], + ] ) + ucs_results = rotor8.run_ucs(synchronous=True) + assert_allclose(ucs_results.wn, exp_rotor_wn) def test_ucs_rotor9(rotor9): - ucs_results = rotor9.run_ucs(num_modes=8) - fig = ucs_results.plot() - expected_x = np.array( + exp_rotor_wn = np.array( + [ + [ + 89.61923784473375, + 120.89655743664413, + 162.59246554800114, + 217.42553100241312, + 287.6459465590863, + 372.935347725518, + 466.2581426686537, + 551.4076269641628, + 613.1237259191993, + 650.4689786203988, + 671.2479800816767, + 682.5353774789179, + 688.6485929814294, + 691.965477524608, + 693.7688181845082, + 694.7506707995042, + 695.2857209934174, + 695.577438207025, + 695.7365318523782, + 695.8233102639163, + ], + [ + 124.8108608009496, + 168.92618197956605, + 228.5753022242499, + 309.1399173060739, + 417.7337051968829, + 563.5536938601297, + 757.9439326311913, + 1013.3462956039914, + 1338.3333635278868, + 1716.48711727428, + 2004.8471556137, + 2078.8761601177744, + 2097.725398329982, + 2104.9250739620065, + 2108.224809182831, + 2109.8706285775806, + 2110.7269388638374, + 2111.182371668304, + 2111.427443127732, + 2111.5601496202053, + ], + [ + 976.610014941276, + 979.8071596437271, + 985.7435440438873, + 996.8720841293019, + 1018.0435385669914, + 1059.040021915284, + 1138.8336540231767, + 1287.1396017447778, + 1530.728131538211, + 1866.6004173519655, + 2224.740955271474, + 2566.1521460703593, + 2812.5568536324695, + 2954.956819348177, + 3030.8828351975917, + 3070.974195701119, + 3092.3284826481, + 3103.806643182663, + 3110.0148616594533, + 3113.3854057718704, + ], + [ + 2159.640747974065, + 2159.756470270719, + 2159.969971035258, + 2160.366035282546, + 2161.1082776995636, + 2162.5260227884364, + 2165.334410630516, + 2171.3131180341056, + 2186.1348008252444, + 2238.1871377763114, + 2490.9613952181917, + 2978.354255657329, + 3456.7805801535656, + 3814.959990456543, + 4040.908191796628, + 4171.396315088025, + 4244.09186723666, + 4284.062958149884, + 4305.937322108033, + 4317.887004191494, + ], + ] + ) + ucs_results = rotor9.run_ucs() + assert_allclose(ucs_results.wn, exp_rotor_wn) + + exp_rotor_wn = np.array( [ - 1.00000000e06, - 1.83298071e06, - 3.35981829e06, - 6.15848211e06, - 1.12883789e07, - 2.06913808e07, - 3.79269019e07, - 6.95192796e07, - 1.27427499e08, - 2.33572147e08, - 4.28133240e08, - 7.84759970e08, - 1.43844989e09, - 2.63665090e09, - 4.83293024e09, - 8.85866790e09, - 1.62377674e10, - 2.97635144e10, - 5.45559478e10, - 1.00000000e11, + [ + 89.61947064, + 120.89741889, + 162.5960395, + 217.44095694, + 287.71229812, + 373.20676063, + 467.22436151, + 554.01174191, + 618.0228653, + 657.38198264, + 679.5050572, + 691.58822264, + 698.15087843, + 701.71684971, + 703.65712249, + 704.71396966, + 705.29001687, + 705.60412456, + 705.77544067, + 705.86888929, + ], + [ + 126.08902147, + 170.65739193, + 230.92107312, + 312.32101392, + 422.05546021, + 569.45142419, + 766.09122048, + 1025.02536231, + 1357.18701994, + 1759.29639866, + 2139.22396755, + 2225.72234144, + 2241.24852987, + 2246.71661805, + 2249.14967666, + 2250.34715925, + 2250.96610863, + 2251.29417729, + 2251.47039346, + 2251.56572169, + ], + [ + 1006.73169597, + 1010.04096437, + 1016.17902002, + 1027.66388689, + 1049.44232444, + 1091.39953214, + 1172.55460038, + 1322.87033841, + 1570.73575288, + 1916.3041552, + 2265.71692139, + 2655.54276614, + 2969.42618151, + 3163.00092959, + 3269.25575144, + 3325.84538219, + 3356.03761948, + 3372.26520948, + 3381.03936768, + 3385.80171304, + ], + [ + 2282.84557629, + 2282.91523796, + 2283.04371506, + 2283.28189556, + 2283.72771552, + 2284.57736115, + 2286.25338975, + 2289.79373785, + 2298.45136454, + 2328.81587039, + 2527.00408181, + 3011.54457196, + 3511.90570401, + 3911.96812784, + 4184.16765863, + 4351.24520897, + 4448.09334563, + 4502.60333569, + 4532.82984364, + 4549.46314886, + ], ] ) - assert_allclose(fig.data[0]["x"], expected_x) + ucs_results = rotor9.run_ucs(synchronous=True) + assert_allclose(ucs_results.wn, exp_rotor_wn, rtol=1e-6) def test_pickle(rotor8):