Damping Coupling Multirotor

[murilloabs]

Description

This pull request implements linear gear-mesh damping in MultiRotor models and refactors the gear coupling representation so that the same geometric coupling matrix is used for both mesh stiffness and mesh damping.

PR #1426 already include damping_ratio in MultiRotor.__init__(). In this PR the default value is swithced to 0 and now this damping_ratio is used for the linear mesh coupling when backlash is disabled.

The mesh-damping coefficient is calculated according to Yi et al. (2019):

$$ c_m(t) = 2\zeta\sqrt{k_m(t)M_{\mathrm{eq}}} $$

where:

• $\zeta$ is the mesh damping ratio;
• $k_m(t)$ is the instantaneous mesh stiffness;
• $M_{\mathrm{eq}}$ is the equivalent mass of the gear pair;
• $c_m(t)$ is the mesh-damping coefficient.

The stiffness and damping contributions are assembled using the same geometric coupling matrix:

$$ K_{\mathrm{mesh}}(t) = k_m(t),\mathtt{coupling_matrix} $$

$$ C_{\mathrm{mesh}}(t) = c_m(t),\mathtt{coupling_matrix} $$

This approach makes the damping coupling analogous to the stiffness coupling and preserves the same interaction pattern between the driving and driven rotors.

Theoretical reference

The mesh-damping formulation follows:

Yi, Y.; Huang, K.; Xiong, Y.; Sang, M. Nonlinear dynamic modelling and analysis for a spur gear system with time-varying pressure angle and gear backlash. Mechanical Systems and Signal Processing, volume 132, pages 18–34, 2019.
DOI: 10.1016/j.ymssp.2019.06.013

The implementation uses the damping expression from this reference and assembles the resulting coefficient into the global damping matrix using the same coupling_matrix used by the mesh stiffness.

The geometric coupling-matrix formulation used by ROSS is based on the existing rotor-dynamics formulation described by Kaplan et al. (2013). The corresponding bibliography entry is:

Kaplan, J.; Dousti, S.; Allaire, P.; Nichols, B.; Dimond, T.; Untaroiu, A. Rotor Dynamic Modeling of Gears and Geared Systems. Proceedings of the ASME Turbo Expo, 2013. DOI: 10.1115/GT2013-94654

A complementary reference for the geared-rotor formulation is Visnadi (2022). This thesis is included here as supporting literature for the use of the gear-coupling formulation in the assembly of mesh-related stiffness and damping contributions. In the present implementation, the same dimensionless geometric coupling_matrix is therefore reused for both terms, while the physical scalar coefficients remain distinct:

$$ K_{\mathrm{mesh}}(t)=k_m(t),\mathtt{coupling_matrix} $$

$$ C_{\mathrm{mesh}}(t)=c_m(t),\mathtt{coupling_matrix} $$

Visnadi, L. B. Efeito de trinca em engrenagens de dente reto na resposta dinâmica do rotor. 2022. Tese (Doutorado em Engenharia Mecânica) — Universidade Estadual de Campinas, Faculdade de Engenharia Mecânica, Campinas, SP, 2022. Disponível em: https://hdl.handle.net/20.500.12733/5500.

New damping-ratio parameter

The MultiRotor class now accepts a damping_ratio parameter:

multi_rotor = MultiRotor(
    rotor1,
    rotor2,
    coupled_nodes=(node_1, node_2),
    gear_mesh_stiffness=gear_mesh_stiffness,
    damping_ratio=0.07,
)

The default value is:

damping_ratio=0.0

This preserves the previous linear numerical behavior because linear mesh damping is disabled when the damping ratio is zero.

The parameter is propagated through the model as follows:

MultiRotor.damping_ratio
        |
        v
Mesh.damping_ratio
        |
        v
Backlash.damping_ratio

The same damping ratio is therefore available to both the linear mesh-damping matrix and the nonlinear backlash-force formulation.

Equivalent mass

The equivalent gear-pair mass is now calculated by Mesh and stored as Mesh.M_eq. It is available whether or not backlash is enabled.

The equivalent mass is calculated as:

$$ M_{\mathrm{eq}} = \frac{I_{p1}I_{p2}} {I_{p2}R_{b1}^{2} + I_{p1}R_{b2}^{2}} $$

where $I_{p1}$ and $I_{p2}$ are the polar moments of inertia and $R_{b1}$ and $R_{b2}$ are the base-circle radii of the two gears.

When backlash is enabled, the same value is passed from Mesh to Backlash:

Backlash(
    ...,
    M_eq=self.M_eq,
    damping_ratio=self.damping_ratio,
)

This avoids duplicating the equivalent-mass calculation and ensures that the linear and nonlinear formulations use the same physical parameter.

Code execution flow

The main initialization flow is:

MultiRotor.__init__()
        |
        | receives damping_ratio
        v
Mesh(...)
        |
        | stores stiffness and damping_ratio
        | calculates M_eq
        v
Mesh.M_eq
        |
        | if backlash is enabled
        v
Backlash(..., M_eq, damping_ratio)
        |
        v
Create coupling_matrix

When a MultiRotor object is created:

1. The driving and driven rotors are received.
2. A Mesh object is created.
3. Mesh stores the mesh stiffness and damping ratio.
4. Mesh.M_eq is calculated.
5. If backlash is enabled, M_eq and damping_ratio are passed to Backlash.
6. The gear degrees of freedom are identified.
7. The geometric coupling_matrix is calculated.
8. The appropriate stiffness and damping coupling functions are selected.

Coupling-matrix refactor

The previous attribute:

MultiRotor.K_coupling

was renamed to:

MultiRotor.coupling_matrix

The matrix itself was not changed. coupling_matrix preserves exactly the same formulation, equations, entries, and numerical structure previously used by K_coupling. This is the existing dimensionless geometric coupling matrix implemented in compute_coupling_matrix(), following the Kaplan et al. (2013) rotor-dynamics formulation already used by ROSS. Only the property name was changed to reflect that the matrix is now used by both stiffness and damping.

Therefore, this pull request does not modify the coupling-matrix formulation. It only reuses the existing matrix for the damping contribution in addition to the stiffness contribution.

The same unchanged matrix defines the diagonal and cross-rotor interaction terms:

Driving-rotor DOFs  <---->  Driven-rotor DOFs
          \              /
           coupling_matrix

Stiffness-matrix flow

When MultiRotor.K() is called, the sequence is:

MultiRotor.K()
        |
        v
Assemble rotor-1 stiffness matrix
        |
        v
Assemble rotor-2 stiffness matrix
        |
        v
Join the rotor matrices
        |
        v
Apply selected mesh-stiffness coupling

For the linear model, the mesh-stiffness contribution is:

$$ K_{\mathrm{mesh}}(t) = k_m(t),\mathtt{coupling_matrix} $$

It is added to the gear degrees of freedom as:

K0[np.ix_(dofs, dofs)] += (
    self.mesh.stiffness * self.coupling_matrix
)

Damping-matrix flow

When MultiRotor.C() is called, the sequence is:

MultiRotor.C()
        |
        v
Assemble rotor-1 damping matrix
        |
        v
Assemble rotor-2 damping matrix
        |
        v
Join the rotor matrices
        |
        v
Calculate c_m
        |
        v
Apply C_mesh()

The coefficient is calculated as:

$$ c_m(t) = 2\zeta\sqrt{k_m(t)M_{\mathrm{eq}}} $$

The damping contribution is then assembled analogously to stiffness:

$$ C_{\mathrm{mesh}}(t) = c_m(t),\mathtt{coupling_matrix} $$

The implementation is conceptually:

c_m = (
    2.0
    * self.mesh.damping_ratio
    * np.sqrt(self.mesh.stiffness * self.mesh.M_eq)
)

C0[np.ix_(dofs, dofs)] += (
    c_m * self.coupling_matrix
)

When damping_ratio=0.0, c_m is zero and no additional cross-rotor damping is added.

Linear and backlash coupling selection

The coupling functions are selected according to the backlash configuration:

if self.mesh.backlash:
    self.add_coupling_stiffness = lambda K0: K0
    self.add_coupling_damping = lambda C0: C0
else:
    self.add_coupling_stiffness = self.K_mesh
    self.add_coupling_damping = self.C_mesh

When backlash is disabled:

K() -> adds K_mesh()
C() -> adds C_mesh()

When backlash is enabled:

K() -> does not add linear K_mesh()
C() -> does not add linear C_mesh()
Backlash -> calculates the nonlinear gear-contact force

This prevents the linear mesh contributions from being added on top of the nonlinear backlash formulation.

Time-varying stiffness and damping

When time-varying mesh stiffness is enabled, Mesh.stiffness is updated before both global matrices are assembled:

Time step t
    |
    v
Calculate k_m(t)
    |
    v
Update Mesh.stiffness
    |
    +--------------------+
    |                    |
    v                    v
Assemble K(t)         Assemble C(t)
    |                    |
    |                    v
    |              Calculate c_m(t)
    |                    |
    |                    v
    |              Add C_mesh(t)
    |
    v
Add K_mesh(t)

The damping coefficient always uses the current mesh stiffness:

$$ c_m(t) = 2\zeta\sqrt{k_m(t)M_{\mathrm{eq}}} $$

Backlash flow

When backlash is active, the gear interaction is calculated through the nonlinear contact-force model:

Relative gear displacement
        |
        v
Evaluate backlash clearance
        |
        +--> Teeth separated:
        |        contact force is zero
        |
        +--> Teeth engaged:
                 calculate nonlinear mesh force
                 include stiffness and damping effects

The damping_ratio is passed to Backlash and used in the nonlinear force formulation. The linear K_mesh and C_mesh matrices are not added separately when backlash is active.

This separates the two formulations:

Linear model:
    gear coupling is represented by K_mesh and C_mesh

Backlash model:
    gear coupling is represented by the nonlinear backlash force

This avoids double counting of mesh stiffness and damping.

Rebuild behavior

The damping_ratio value is preserved whenever the MultiRotor model is rebuilt. This applies to operations such as adding nodes, adding elements, and reconstructing the internal rotor model.

Main code changes

• Added damping_ratio to MultiRotor.
• Propagated damping_ratio from MultiRotor to Mesh.
• Propagated damping_ratio from Mesh to Backlash.
• Added Mesh.M_eq, available with or without backlash.
• Moved the equivalent-mass calculation from Backlash to Mesh.
• Renamed K_coupling to coupling_matrix.
• Added the C_mesh() method.
• Updated K_mesh() to use coupling_matrix.
• Updated C() and K() to apply the selected coupling behavior.
• Updated the time-varying stiffness integration path.
• Preserved damping_ratio during rebuild operations.
• Added tests for equivalent mass, linear damping, cross-rotor terms, time-varying stiffness, backlash behavior, and rebuilds.
• Updated the MultiRotor, Mesh, Backlash, and C_mesh() docstrings.
• Updated the multirotor tutorial.
• Updated the existing bibliography entry for Yi et al. (2019).
• Did not modify the graphical interface or its Python exporter.

Compatibility notes

This change introduces a breaking API rename:

MultiRotor.K_coupling

was replaced by:

MultiRotor.coupling_matrix

No compatibility alias was added.

The default value is:

damping_ratio=0.0

This disables linear mesh damping and preserves the previous linear numerical behavior.

[codecov-commenter]

⚠️ Please install the to ensure uploads and comments are reliably processed by Codecov.

Codecov Report

❌ Patch coverage is 95.45455% with 1 line in your changes missing coverage. Please review.
✅ Project coverage is 82.81%. Comparing base (4aa371d) to head (bed5842).
⚠️ Report is 1 commits behind head on main.

Files with missing lines	Patch %	Lines
ross/multi_rotor/multi_rotor.py	93.75%	1 Missing ⚠️
❗ Your organization needs to install the Codecov GitHub app to enable full functionality.

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1427      +/-   ##
==========================================
+ Coverage   82.59%   82.81%   +0.22%     
==========================================
  Files         173      173              
  Lines       27549    27561      +12     
==========================================
+ Hits        22754    22826      +72     
+ Misses       4795     4735      -60     

Files with missing lines	Coverage Δ
ross/multi_rotor/mesh.py	63.75% <100.00%> (+0.11%)	⬆️
ross/units.py	90.41% <ø> (ø)
ross/multi_rotor/multi_rotor.py	89.58% <93.75%> (+1.45%)	⬆️

... and 1 file with indirect coverage changes

---

Continue to review full report in Codecov by Harness.

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

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.