Incorporate time integration to account for speed variations

[jguarato]

This issue was created in order to document and detail the task of implementing the numerical integrator that was assigned to me. This task may resolve other issues: #429, #884, #1001.

---

Currently, it is not possible to solve the transient model with the rotor speed varying over time in ROSS, nor to incorporate other information into the model that depends on time. One solution would be to use a numerical integrator to solve the following equation of motion for a rotor:

$$ M\ \ddot{q}+[C+ G\ \Omega(t)]\ \dot{q} + [K + K_{st}\ \dot{\Omega}(t)]\ {q} = F $$

The function solve_ivp() available in SciPy integrate library can be applied in this case. This function numerically integrates a set of first-order vector ordinary differential equations (ODEs):

$$ \frac{d\mathbf{y}}{dt} = \mathrm{f}(t,\ \mathbf{y}) $$

given the initial conditions $\mathbf{y}(0)=y_0$, and it also offers 6 different integration methods. Then, it would be possible to choose a more appropriate method for each type of problem, allowing to solve, for example, non-stiff problems using explicit methods (“RK23”, “RK45”, “DOP385”), and also stiff problems using implicit methods (“Radau”, “BDF”).

A higher-order ODE can be reduced to a first-order differential equation by introducing intermediate derivatives into the $\mathbf{y}$ vector. So, this implies that:

$$ \mathbf{y} = [q, \dot{q}] $$

and accordingly $\mathrm{f}(t,\ \mathbf{y})$ must return:

$$ \frac{d\mathbf{y}}{dt} = [\dot{q}, \ddot{q}] $$

where the second-order derivative could be calculated from the equation of motion by inverting the $M$ matrix as:

$$ \ddot{q} = M^{-1} [F - [C+ G\ \Omega(t)]\ \dot{q} - [K + K_{st}\ \dot{\Omega}(t)]\ {q} ] $$

The values of rotational speed $\Omega(t)$ and acceleration $\dot{\Omega}(t)$ varying over time can be stored in arrays. Such arrays must be compatible with the array of times in which the computed solution will be stored.

Consequently, the integration can be done by calling:

scipy.integrate.solve_ivp(fun, t_span, y0, method, t_eval, args)

and passing the callable function that will be integrated fun(t, y), the interval of integration t_span=[t0, tf] and the initial state y0. It is optional, but it could be important to pass the integration method, the vector of times in which the computed solution will be stored t_eval, and a tuple of additional arguments that could be passed to fun(t, y, *args).

At first, I defined the function fun(t, y, *args) to be integrated as follows:

def _equation_of_motion(self, t, y, t0, dt, speed, accel, inv_M, C, G, K, Kst, F):
	
  size = int(len(y) / 2)
  disp_resp = y[:size]
  velc_resp = y[size:]

  i = round((t - t0) / dt)
  C_aux = C + G * speed[i]
  K_aux = K + Kst * accel[i]
  F_total = F[:, i]

  accl_resp = inv_M @ (F_total - (C_aux @ velc_resp) - (K_aux @ disp_resp))

  dydt = np.zeros(2 * size)
  dydt[:size] = velc_resp
  dydt[size:] = accl_resp

  return dydt

Here, I would like to highlight two points:

• Accessing the rotational speed speed[i] at a specific time t can be done by a conversion based on the initial time t0 and time step dt to obtain the respective array index.
• Passing the model matrices as arguments makes the function generic, so it could be used to solve the model in both modal and physical spaces. However, this may change later if needed.

When choosing to solve in the modal space, there is a possibility to reduce model based on the selected number of modes. A reduced model implies a considerable reduction in computational cost. Another option, suggested by Prof. Aldemir, to obtain the reduced model would be to pass the frequency range of interest.

[jguarato]

I tested a simple case to compare the numerical integration considering the representation of the model in both spaces. For the modal space, just the first 12 modes were considered in order to obtain the reduced model. For both spaces, the “RK45” integration method was considered. This test case is the same example described in run_time_response():

# Rotor without damping with 6 shaft elements, 2 disks and 2 bearings
rotor = rotor_example()
speed = 500.0
size = 1000
node = 3
probe1 = (3, 0)
t = np.linspace(0, 10, size)
F = np.zeros((size, rotor.ndof))
F[:, 4 * node] = 10 * np.cos(2 * t)
F[:, 4 * node + 1] = 10 * np.sin(2 * t)
response = rotor.run_time_response(speed, F, t)
fig1 = response.plot_1d(probe=[probe1])

Running the code with time integration in the physical space took a total of 129 s, while in the modal space it took a total of 5.71 s. Therefore, considering the reduced model it was possible to reduce ~96% of the execution time. The responses for node 3 are compared and presented below:

The response obtained by integrating in modal space seems to be very close to the response obtained by lsim().

In terms of computational cost, I think it would be better to implement the integration solver in modal space. The idea would be to use a single solver for Rotor and Defect classes.

@raphaeltimbo, any suggestions for additional testing? Which direction should I take with the implementations?

[raphaeltimbo]

Hi @jguarato !
I was wondering why the 'physical space' results are showing higher values of vibration at the end of the plot. Maybe results are diverging?
If you could please increase the time array so that we could evaluate that.
Also, if you could provide time results for the run_time_response as well.

[jguarato]

@raphaeltimbo, you are right. I increased the time array and the results for the ‘physical space’ diverged using the RK45 integration method.

I found the following article that may explain what happened in this case: Stability of Numerical Integration Techniques for Transient Rotor Dynamics. The author analyzed the numerical stability of some integration techniques, including the Runge-Kutta method, applied to transient rotor dynamics. One of the author's conclusions was that the high frequency mode usually determines the maximum time step for a stable solution. Therefore, considering an undamped system, the Runge-Kutta method is unstable for a time step greater than approximately:

$$ \Delta t \gtrapprox \frac{3}{\omega_n} $$

where $\omega_n$ is the highest natural frequency.

Consequently, calculating for the problem at hand, this method is unstable for $\Delta t \gtrapprox 6 \times 10 ^{-5}$ s. However, in the previously tests, I considered $\Delta t = 10 ^{-3}$ s. Thus, reducing the time step to $\Delta t = 10 ^{-5}$ s, the obtained results were:

In this sense, an implicit method such as BDF ends up being more suitable for this case, for which it is possible to use a larger time step to obtain a stable solution:

We can then compare these solutions with the results of the run_time_response() function:

In the case of solving the equations in modal space, as the reduced model was considered based on the first 12 modes, it was possible to obtain a stable solution for a larger time step. Therefore, solving the equations in modal space with RK45 and without reducing the model also leads to an unstable solution, as can be seen below:

[raphaeltimbo]

To implement the newmark integration:

• Add a newmark function to the utils.py module;
• Import the newmark function in the rotor_assembly.py with from ross.utils import newmark;
• In the .run_time_response method, add an argument integrator where the user can set integrator='newmark'.

[ross-bott]

Hi there!
I have marked this issue as stale because it has not had activity for 45 days.
Consider the following options:

• If the issue refers to a large task, break it in smaller issues that can be solved in
  less than 45 days;
• Label the issue as wontfix or wontfix for now and close it.

[jguarato]

Discussions ensued regarding the previously acquired results, leading to the decision to adopt the Newmark method over the Scipy integrator. The new implementations have been incorporated into the PR #1032. The results of the tests conducted to validate the implementations can also be found in the comments of the PR.