Add Newmark integrator in time response method

[jguarato]

Resolves #1021

---

Hi @raphaeltimbo,

I created this pull request to report what I have been doing regarding the numerical integrator.
If there's anything I need to change, feel free to tell me.
I have some questions related to some of the implementations, and I will describe them better in the comments below.

[jguarato]

Description

The Newmark integrator was implemented to take into account the speed variation in the rotor time response method. It can be used by passing an argument, integrator="newmark", to run_time_response() or passing speed as an array of values with the same length as t (time). The source of the new integrator is now available in utils.py (4 new functions were created).

The pseudo-modal method was also implemented to be considered in numerical integration, and it is only activated if the num_modes (number of modes) argument is passed to run_time_response(). In this method, the model is solved in the modal domain and can be reduced according to the first $n$ chosen modes. If num_modes is left out, the direct method is considered and the model is solved in the physical domain without reduction.

Doubts

I am not sure if what I did on lines 989-1011 of utils.py is correct. To reduce the computational time, I tried to separate the bearing matrices from the other elements, and update them with the respective speed of the integration time step. This implementation was made for the cases where the bearing coefficients vary with speed, and the speed varies over time.

What do you think?

Tests

Some tests were carried out to evaluate the integrator response and are presented below.

[jguarato]

Test 1

Model (with proportional damping):
• 36 elements
• 2 disks
• 2 bearings (coefficients do not vary with speed and there are no cross coefficients)

Analysis in the transient regime:
• Max speed = 840 rpm (87.96 rad/s)
• Unbalance at nodes [12, 23]
 ◦ Magnitude = [0.00034, 0.00034] kg.m
 ◦ Phase = [70, 30] degrees
• Probes at nodes [3, 34]

After assembling the rotor, the simulation was run by making the following call:

response = rotor.run_time_response(speed, F, t, integrator="newmark")

The results were compared against the results obtained with the code of LMEst - UFU:

Then, the reduced model was assumed by the pseudo-modal method for the first 12 modes:

response = rotor.run_time_response(speed, F, t, integrator="newmark", num_modes=12)

The results were compared against the results obtained with the direct method:

In this case, the pseudo-modal method seemed to present results very similar to the direct method. The time taken to run both simulations was 243 s for the direct method (physical domain) and 22 s for the pseudo-modal method with 12 modes (approximately 91% reduction in computational cost).

[jguarato]

Test 2

Model:
• 30 elements
• 13 disks
• 2 bearings (coefficients vary with speed)

Analysis in the transient regime:
• Max speed = 1800 rpm (188.5 rad/s)
• Unbalance at nodes [8, 16]
 ◦ Magnitude = [0.0832, 0.1097] kg.m
 ◦ Phase = [75, 330] degrees
• Probes at nodes [3, 21]

For this case, the results were obtained by reducing the model to the first 12, 48, and 72 modes, which were then compared with the direct method:

As seen from the graphs above, considering only the first 12 modes in the pseudo-modal method is not representative in this case, as evidenced by the results of test 1. One possible explanation is the significant magnitude of the bearing cross coefficients. Therefore, to improve results, a larger number of modes must be considered, but this leads to an increase in computational time. For instance, the direct method required 3265 s to execute, whereas the pseudo-modal method took 1288 s with 12 modes, 2803 s with 48 modes, and 3125 s with 72 modes.

---

⏱️ Update

I finished the comparisons with the LMEst code, and with XLTRC:

For XLTRC I ran 3 different maximum speeds:

[codecov-commenter]

Codecov Report

Attention: Patch coverage is 79.33884% with 25 lines in your changes are missing coverage. Please review.

Project coverage is 85.28%. Comparing base (4b78ef6) to head (ddf46d1).

❗ Your organization needs to install the Codecov GitHub app to enable full functionality.

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1032      +/-   ##
==========================================
- Coverage   85.32%   85.28%   -0.05%     
==========================================
  Files          35       35              
  Lines        7620     7735     +115     
==========================================
+ Hits         6502     6597      +95     
- Misses       1118     1138      +20     

Files	Coverage Δ
ross/rotor_assembly.py	93.38% <80.00%> (+0.27%)	⬆️
ross/utils.py	91.46% <79.27%> (-4.36%)	⬇️

---

Continue to review full report in Codecov by Sentry.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 4b78ef6...ddf46d1. Read the comment docs.

[jguarato]

Additional

Resolves #884

I added a new argument that can be passed to the integrator add_to_RHS:

"""
add_to_RHS : callable, optional
    An optional function that computes and returns an additional array to be added to
    the right-hand side of the equation of motion. This function should take arguments
    corresponding to the current state of the rotor system, including the time step
    number, displacements, and velocities. It should return an array of the same length
    as the degrees of freedom of the rotor system (`rotor.ndof`). This function allows
    for the incorporation of supplementary terms or external effects in the rotor system
    dynamics beyond the specified force input during the time integration process.
"""

def my_callable_function(step, displacement, velocity):

    ...

    return ext_force

	
response = rotor.run_time_response(speed, F, t, integrator="newmark", add_to_RHS=my_callable_function)