# Dynamic substructuring based on the complex Udwadia-Kalaba formalism

François Fabre, Jean-Loïc Le Carrou, Baptiste Chomette

► To cite this version:

François Fabre, Jean-Loïc Le Carrou, Baptiste Chomette. Dynamic substructuring based on the complex Udwadia-Kalaba formalism. International Conference on Noise and Vibration Engineering, Sep 2024, Louvain (Leuven), Belgium. hal-04942911

HAL Id: hal-04942911

https://hal.science/hal-04942911v1

Submitted on 12 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Dynamic substructuring based on the complex Udwadia-Kalaba formalism

F. Fabre*<fnref n="1" />, J.-L. Le Carrou<fnref n="1" />, B. Chomette<fnref n="2" />

## Abstract

Modal dynamic substructuring is widely used to perform structural dynamic analysis. Recently, an efficient modal formulation using (free interface) real modes of the substructures has been proposed based on the Udwadia-Kalaba formalism. The latter expresses the subsystems' constrained response from their unconstrained response by explicitly solving the constraining forces. From an experimental point of view, in presence of non-classical damping, the modal formulation based on real modes presents certain difficulties due to its use of the modal damping matrix. The study presents an alternative formulation based on complex modes of the subsystems, requesting only the knowledge of vibration poles and complex mode shapes. As an example of application, the dynamic response of three harps from Central Africa exhibiting strongly complex modes is simulated by assuming both a proportional and non-proportional damping. Significant differences of frequency content are observed due to the absence of coupling by damping between some sets of modes.

## 1 Introduction

When analyzing the dynamical behavior of a large and/or complex system, it may be convenient to consider it as a multibody structure and performing the analysis at the substructure level. The behavior of the global system is then obtained via coupling between the substructures, described by constraint conditions. This approach, known as dynamic substructuring, is also particularly attractive when dealing with hybrid assemblies consisting in both theoretically and experimentally characterized subsystems.

The field of dynamic substructuring was first introduced to deal with systems made of rigid components. During the last 60 years, methods allowing to account for the visco-elasticity of the substructures have been developed. A classification of the existing methods was proposed by de Klerk et al. [1] who grouped them in three categories: the physical, modal and frequency domains.

A review of some famous physical domain methods was proposed by Lalausa and Bachau [2] including for example Maggi's method [3, 4] and the Udwadia-Kalaba (U-K) formulation [5, 6, 7, 8, 9]. Physical domain formulations require the knowledge of physical mass, damping and stiffness matrices of each subsystem, making them unpractical in the frame of experimental substructuring.

On the other hand, frequency based substructuring (FBS) formulations, such as the admittance [10, 11], impedance [12] and the Lagrange multipliers FBS [13, 14] coupling, only require the knowledge of frequency response functions, whose measurement present no particular difficulty but implies that the substructures are linear time invariant (LTI) in steady-state.

Finally, modal domain formulations, also known as component mode synthesis (CMS), are based on a modal reduction of the substructures which makes them not as straight forward to use as FBS but offers a solution

<footnote n="1">¹ Sorbonne Université, CNRS, Institut Jean Le Rond d'Alembert, Équipe Lutheries-Acoustique-Musique, F-75005 Paris, France</footnote>

<footnote n="2">² Ecole Centrale de Lyon, CNRS, ENTPE, LTDS, UMR5513, 69130 Ecully, France</footnote>

<footnote>*now working at PIMM, CNRS, Arts et Métiers, Institut Polytechnique de Paris, Palaiseau, France</footnote>
for analyzing non-LTI multibody systems. An overview of CMS methods based on all kinds of component modes was proposed by Craig [15]. Examples of such methods are the Craig-Brampton [16] (combination of constraint modes and fixed-interface normal modes), MacNeal [17] and Rubin [18] (combination of attachment modes and free-interface normal modes) methods. In the last decade, in the frame of musical acoustics, Antunes and Debut [19] adapted the U-K formalism to continuous flexible systems using free-interface normal modes. Its relevance has since been demonstrated when dealing with non-linearities [20, 21, 22, 23] and intermittent contacts [24]. Their modal formulation allows for redundant and non-ideal constraints due to the U-K formalism and circumvents the need for knowledge of physical mass, damping and stiffness matrices of the substructures.

While the modal U-K formulation proposed by Antunes and Debut makes using the U-K formalism practical for experimental substructuring, certain limitations persist when dealing with non-classically damped substructures. Indeed, the decomposition of the substructures' response on their normal basis requires to know their modal damping matrix whose estimation is still nowadays a difficult task when this matrix is non-diagonal. To remedy this, this paper presents an alternative modal U-K formulation based on complex modes of the dissipative substructures.

The original U-K formulation in the physical domain is recalled in section 2.1, followed by the presentation of an alternative formulation in the modal domain based on the decomposition of each constrained substructure response on its unconstrained (free interface) complex modal basis. The advantages of this alternative formulation are then illustrated in section 2.2 using experimental data collected on a corpus of harps from Central Africa.

# 2 Dynamic substructuring using the Udwadia-Kalaba formalism

## 2.1 Physical degrees of freedom

Let $y(t)$ represent the response, and $x(t)$ the degrees of freedom, of a reciprocal discrete multibody mechanical system which consists of $J$ subsystems constrained via constraining forces $f_c$ and subjected to external constraint-independent forces $f_e$. Assuming that each subsystem has its degrees of freedom $x^j$ verifying a second order model of the form

$$ \begin{gathered} \mathbf{M}^j \ddot{\mathbf{x}}^j(t) + \mathbf{C}^j \dot{\mathbf{x}}^j(t) + \mathbf{K}^j \mathbf{x}^j(t) = (\mathbf{D}^j)^T (\mathbf{f}_e^j(\mathbf{x}, \dot{\mathbf{x}}, t) + \mathbf{f}_c^j(\mathbf{x}, \dot{\mathbf{x}}, t)) \\ \mathbf{y}^j(t) = \mathbf{D}^j \mathbf{x}^j(t) \end{gathered} \quad (1) $$

with $\mathbf{M}^j$, $\mathbf{C}^j$ and $\mathbf{K}^j$ its mass, damping and stiffness matrices, $\mathbf{D}^j$ its input/output shape matrix, $\mathbf{f}_e^j$ external constraint-independent forces and $\mathbf{f}_c^j$ constraining forces.

Assembling Eq. (1) for all subsystems into a global system leads to

$$ \begin{gathered} \mathbf{M}\ddot{\mathbf{x}}(t) + \mathbf{C}\dot{\mathbf{x}}(t) + \mathbf{K}\mathbf{x}(t) = \mathbf{D}^T (\mathbf{f}_e(\mathbf{x}, \dot{\mathbf{x}}, t) + \mathbf{f}_c(\mathbf{x}, \dot{\mathbf{x}}, t)) \\ \mathbf{y}(t) = \mathbf{D}\mathbf{x}(t) \end{gathered} \quad (2) $$

where the assembled matrices and vectors of Eq. (2) are defined

$$ \begin{gathered} \mathbf{x} = \begin{Bmatrix} \mathbf{x}^1 \\ \vdots \\ \mathbf{x}^J \end{Bmatrix}, \quad \mathbf{f}_e = \begin{Bmatrix} \mathbf{f}_e^1 \\ \vdots \\ \mathbf{f}_e^J \end{Bmatrix}, \quad \mathbf{f}_c = \begin{Bmatrix} \mathbf{f}_c^1 \\ \vdots \\ \mathbf{f}_c^J \end{Bmatrix}, \quad \mathbf{M} = \operatorname{diag}(\mathbf{M}^1, \dots, \mathbf{M}^J) = \begin{bmatrix} \mathbf{M}^1 & & & \mathbf{0} \\ & \ddots & & \\ \mathbf{0} & & & \mathbf{M}^J \end{bmatrix}, \\ \mathbf{C} = \operatorname{diag}(\mathbf{C}^1, \dots, \mathbf{C}^J), \quad \mathbf{K} = \operatorname{diag}(\mathbf{K}^1, \dots, \mathbf{K}^J), \quad \mathbf{D} = \operatorname{diag}(\mathbf{D}^1, \dots, \mathbf{D}^J). \end{gathered} $$

Constraining forces can be expressed through Lagrange multipliers [25] $\lambda$ as

$$ \mathbf{f}_c(\dot{\mathbf{x}}, \mathbf{x}, t) = -\mathbf{A}(\dot{\mathbf{x}}, \mathbf{x}, t)^T \lambda(t) \quad (3) $$
where $A(\dot{x}, x, t)$ is the constraint matrix, associated with $b(\dot{x}, x, t)$ a vector function of the motion, corresponding to the following system of P holonomic and non-holonomic constraints in terms of accelerations

$$A(\dot{x}, x, t) \ddot{y} = b(\dot{x}, x, t). \quad (4)$$

For the sake of readability, the dependencies on $\dot{x}$, $x$ and $t$ are omitted in the following.

Combining Eqs. (2–4) leads to the following augmented differential-algebraic equation (DAE) system

$$\begin{bmatrix} M & (AD)^T \\ AD & 0 \end{bmatrix} \begin{Bmatrix} \ddot{x} \\ \lambda \end{Bmatrix} = \begin{Bmatrix} M\ddot{x}_u \\ b \end{Bmatrix} \\ y = Dx \quad (5)$$

where $\ddot{x}_u = M^{-1}(D^T f_e - C\dot{x} - Kx)$ represents the unconstrained response of the subsystems.

Solving Eq. (5), in the least square sense [26] since $A$ may be rank deficient, gives the differential equation verified by the dynamic response of the constrained system

$$\begin{gathered} \ddot{x} = \ddot{x}_u + M^{-1}\Delta(b - \tilde{A}\ddot{x}_u) \\ y = Dx \end{gathered} \quad (6)$$

with $\tilde{A} = AD$ and $\Delta = \tilde{A}^T(\tilde{A}M^{-1}\tilde{A}^T)^\dagger$.

Finally, Eq. (6) may be further simplified by noticing that $\Delta = M^{1/2}(\tilde{A}M^{-1/2})^\dagger = M^{1/2}B^\dagger$ ($\dagger$ denoting the Moore-Penrose pseudo-inverse), so that the main result from Udwadia and Kalaba [5, 7] is

$$\begin{gathered} \ddot{x} = \ddot{x}_u + M^{-1/2}B^\dagger(b - \tilde{A}\ddot{x}_u) \\ y = Dx. \end{gathered} \quad (7)$$

The elegance of Eq. (7) lies in the expression of the constrained response of the discrete multibody system in terms of the unconstrained responses of its subsystems.

From a practical point of view, Eq. (7) has the disadvantage of requesting the knowledge of the physical mass, damping and stiffness matrices whose estimations are not trivial when dealing with experimental data. Moreover, these matrices may be fully populated or at least far from sparse so that solving Eq. (7) is not computationally efficient. To remedy this, Antunes & Debut [19] recently adapted the U-K formulation by means of a modal expansion on the modal bases of associated conservative subsystems and obtained a differential equation equivalent to Eq. (7) whose coefficients involve modal mass, damping and stiffness matrices. Their formulation is particularly attractive in the frame of experimental substructuring when the subsystems yield a classical damping since all modal matrices become diagonal and can be relatively easily estimated via experimental modal analysis. However, if the $J$ subsystems do not verify the assumption of classical damping, the modal damping matrix is non-diagonal and its experimental estimation is difficult. The aim of the next section is to present an alternative modal Udwadia-Kalaba formulation allowing to circumvent the limitations induced by the modal damping matrix requirement.

## 2.2 Modal degrees of freedom

In this section, a modal formulation based on complex modes of the subsystems is presented in order to obtain a diagonal system of unconstrained modal coordinates even for non-classically damped substructures.

Starting from the usual modal expansion on the basis of modal coordinates of the dissipative system is written

$$\begin{Bmatrix} x \\ \dot{x} \end{Bmatrix} = \Upsilon q \quad (8)$$
where, for $j = 1, ..., J$ constrained subsystems, the vectors that assemble the corresponding modal responses, as well as the matrices that assemble the eigenvalues and the eigenvectors, are defined

$$ 
\begin{aligned}
\mathbf{q} &\equiv \left\{ \begin{matrix} q^1 \\ \vdots \\ q^J \end{matrix} \right\}, \quad \mathbf{\Psi} \equiv \mathrm{diag}(\mathbf{\Psi}^1, \dots, \mathbf{\Psi}^J), \quad \mathbf{\Lambda} \equiv \mathrm{diag}(\mathbf{\Lambda}^1, \dots, \mathbf{\Lambda}^J), \\
\mathbf{q} &\equiv \left\{ \frac{\mathbf{q}}{\|\mathbf{q}\|} \right\}, \quad \mathbf{\Psi} \equiv [\mathbf{\Psi} \quad \mathbf{\Psi}^\dagger], \quad \mathbf{\Lambda} \equiv \begin{bmatrix} \mathbf{\Lambda} & 0 \\ 0 & \mathbf{\Lambda}^\dagger \end{bmatrix}, \quad \mathbf{\Upsilon} \equiv \begin{bmatrix} \mathbf{\Psi} \\ \mathbf{\Psi} \mathbf{\Lambda} \end{bmatrix}
\end{aligned} 
\quad (9) $$

with, for each subsystem, the modal basis of $2N^j$ unconstrained modes containing $\mathbf{q}^j$, $\mathbf{\Psi}^j$, $\mathbf{\Lambda}^j$ (and their complex conjugate), respectively, the vector of modal coordinates, the complex mode shapes matrix and the matrix of associated poles. The total number of pairs of complex conjugate modes is thus $N_{\text{tot}} = \sum_{j=1}^{J} N^j$.

Eq. (6) may be augmented into the following state-space form

$$ 
\begin{aligned}
& \underbrace{\begin{bmatrix} \mathbf{C} & \mathbf{M} \\ \mathbf{M} & \mathbf{0} \end{bmatrix}}_{\mathbb{U}} \left\{ \begin{matrix} \dot{\mathbf{x}} \\ \ddot{\mathbf{x}} \end{matrix} \right\} - \underbrace{\left( \begin{bmatrix} -\mathbf{K} & \mathbf{0} \\ \mathbf{0} & \mathbf{M} \end{bmatrix} + \begin{bmatrix} \mathbf{ZK} & \mathbf{ZC} \\ \mathbf{0} & \mathbf{0} \end{bmatrix} \right)}_{\mathbb{A}} \left\{ \begin{matrix} \mathbf{x} \\ \dot{\mathbf{x}} \end{matrix} \right\} = \left\{ (\mathbb{I}_{N_{\text{dof}}} - \mathbf{Z}) \mathbf{D}^T \mathbf{f}_e + \mathbf{\Delta} \mathbf{b} \right\} \\
& \mathbf{y} = \begin{bmatrix} \mathbf{D} & \mathbf{0} \end{bmatrix} \left\{ \begin{matrix} \mathbf{x} \\ \dot{\mathbf{x}} \end{matrix} \right\}
\end{aligned} 
\quad (10) $$

with $\mathbf{Z} = \mathbf{\Delta}\tilde{\mathbf{A}}\mathbf{M}^{-1}$ and $\mathbb{I}_{N_{\text{dof}}}$ the identity matrix of size $N_{\text{dof}}$, the total number of dofs of the subsystems.

Then projecting Eq. (10) on the assembled modal basis of the unconstrained dissipative subsystems by pre-multiplying by $\mathbf{\Upsilon}^T$ leads to

$$ \mathbf{\Pi} \dot{\mathbf{q}} - (\mathbf{\Pi} \mathbf{\Lambda} - \mathbf{\Psi}^T \mathbf{\Delta} \tilde{\mathbf{A}} \mathbf{\Psi} \mathbf{\Lambda}^2) \mathbf{q} = \mathbf{\Psi}^T (\mathbb{I}_{N_{\text{dof}}} - \mathbf{\Delta} \tilde{\mathbf{A}} \mathbf{M}^{-1}) \mathbf{D}^T \mathbf{f}_e + \mathbf{\Psi}^T \mathbf{\Delta} \mathbf{b} \quad (11) $$

using the orthogonality properties [27] $\mathbf{\Upsilon}^T \mathbb{U} \mathbf{\Upsilon} = \mathbf{\Pi}$ and $\mathbf{\Upsilon}^T \mathbb{A} \mathbf{\Upsilon} = \mathbf{\Pi} \mathbf{\Lambda}$ (with $\mathbf{\Pi}$ a diagonal normalization matrix equivalent to the modal mass matrix for principal coordinates).

In addition, using the former orthogonality properties, it can be shown that

$$ \mathbf{\Psi}^T \mathbf{\Delta} = \mathbf{\Psi}^T \tilde{\mathbf{A}}^T (\tilde{\mathbf{A}} \mathbf{\Psi} \mathbf{\Pi}^{-1} \mathbf{\Lambda} \mathbf{\Psi}^T \tilde{\mathbf{A}}^T)^\dagger. \quad (12) $$

Defining now the modal constraint matrix $\mathbf{A} = [\mathbf{A} \quad \mathbf{\Lambda}\mathbf{A}] = \tilde{\mathbf{A}}\mathbf{\Psi}$ and $\mathbf{B} = \mathbf{A} (\mathbf{\Pi}^{-1}\mathbf{\Lambda})^{1/2}$, Eq. (12) may be rewritten

$$ \mathbf{\Psi}^T \mathbf{\Delta} = (\mathbf{\Pi}^{-1} \mathbf{\Lambda})^{-1/2} \mathbf{B}^T (\mathbf{B}\mathbf{B}^T)^\dagger. \quad (13) $$

Finally, combining Eqs. (11) and (13) and left multiplying by $\mathbf{\Pi}^{-1}$, a modal formulation of the Udwadia-Kalaba formalism (equivalent to Eq. (7)) based on the unconstrained (free interface) complex modes of the subsystems is found

$$ 
\begin{aligned}
\dot{\mathbf{q}} &= \dot{\mathbf{q}}_u + (\mathbf{\Pi}\mathbf{\Lambda})^{-1/2} \mathbf{B}^T (\mathbf{B}\mathbf{B}^T)^\dagger (\mathbf{b} - \mathbf{A}\mathbf{\Lambda}\dot{\mathbf{q}}_u) \\
\mathbf{y} &= \begin{bmatrix} \mathbf{D} & \mathbf{0} \end{bmatrix} \mathbf{\Upsilon} \mathbf{q}.
\end{aligned} 
\quad (14) $$

where the unconstrained complex modal response is $\dot{\mathbf{q}}_u = \mathbf{\Lambda}\mathbf{q} + \mathbf{\Pi}^{-1}\mathbf{\Psi}^T\mathbf{D}^T\mathbf{f}_e$.

In summary, Eq. (14) allows to compute the constrained dynamic response of reciprocal multibody mechanical systems with non-classical damping and positive definite mass and stiffness matrices (at the subsystems level).
## 3 Application to an experimental example: a corpus of harps from Central Africa

In this section, the advantage of the complex modal U-K formulation is demonstrated by studying a corpus of harps from Central Africa (see Figure 1a) whose experimental modal analysis has shown the presence of strongly complex modes as will be discussed in the following.

The experimental setup used to measure the frequency response functions is shown in Figure 1b. The sensors consists in three accelerometers (PCB M352C65): two at the same position on the neck but oriented into two transverse polarizations and one mounted on the tailpiece measuring its transverse displacement. An automatic hammer (force sensor PCB 086E80) impacts the harps at locations limited to the string attachments points on the tailpiece and the neck in order to limit the duration of the measurements and, doing so, reducing the risks of changes in room temperature and air humidity during the measurements, as those have been observed to induce shifts in modal frequencies. To help further reduce the changes in ambient conditions, an air humidifier is regularly filled with water and turned on during the measurements. The latter do not ensure a constant air humidity in the room but allows to sufficiently slow down the process during the two to three hours of measurement.

<figure><img src="image_2.png" /><figcaption>Figure 1: (a) Corpus of harps included in the measurement. (b) Experimental setup used to measure acceleration responses at string attachment points. Ⓐ: accelerometers, Ⓑ: automatic impact hammer, Ⓒ: bungee rope, Ⓓ: strip of felt.</figcaption></figure>

For each harp, the potential of coupling between the estimated poles is then quantified using the separation criterion proposed by Hasselman [28]

$$ \delta_{jk} = \frac{2\xi_j\omega_j}{|\omega_j - \omega_k|} \ll 1 \qquad (15) $$

with $\xi_j$ the loss factor of mode $j$, $\omega_j$ and $\omega_k$ the undamped natural frequencies of modes $j$ and $k$.

This criterion signifies that the coupling due to non-proportional damping can be neglected if the cross-modal impedance is high. As suggested by Balmès [27], this criterion can be used to identify groups of modes verifying the assumption of block proportional damping and then, for each group, a proper basis is approximated by resolution of an algebraic Riccati equation [29].

The results are shown in Figure 2 for the frequency range [0 - 1000] Hz, modes are sorted by increasing modal frequencies. To help the comparison between the harps, groups of modes which can not be assumed uncoupled by damping (according to a threshold of 0.3 for the separation criterion $\delta_{jk}$) are indicated by blue dashed squares.
<figure><img src="image_3.png" /><figcaption>Figure 2: Separation criterion for the experimentally identified modes of harps of the corpus (frequency range: [0 1000] Hz). Blue dashed squares indicate complete sets of modes (negligible inter-sets coupling) for a threshold of 0.3 for $\delta_{jk}$ (see Eq. (15)). (a) Harp Fg; (b) Harp Fg_2; (c) Harp Nko.</figcaption></figure>

The comparison of the separation criterion between the three harps allows to draw some conclusions. First, on all harps, there appears to be modes strongly coupled by damping with a separation criterion higher than 0.7 as can be seen, for example, on Figure 2b between modes 20 to 24 and 25 to 30. This indicates that the vibration of harp bodies will be erroneously modeled by a modal model with proportional damping and amplitude differences between the modal Udwadia-Kalaba based on real modes and complex modes should thus be expected for all harps.

In order to test this hypothesis, physical simulations are conducted, for each harp, assuming both a proportionally and non-proportionally damped model. These simulations are all obtained using the complex modal U-K formulation presented in section 2.2. The simulated signal consists in the acceleration response of the soundboard, at the attachment point of the 5th string (tuned at 247 Hz), to a pluck applied at 3/5th of the string length and modeled as a ramp of force of maximum amplitude 1.1 N and duration 1 s. The sampling frequency is set at 800 kHz. The string parameters, shown in Table 1, are kept for all harps so that only modal parameters of the harp bodies vary. The string termination is imposed as perfectly rigid at the end opposite the string/tailpiece coupling point by choosing its unconstrained (free interface) modal basis as that of a fixed-free ideal string. Its motion is described using 150 modes whose modal damping is experimentally estimated. The simulated waveforms and spectra for the three harps are shown in Figure 3.

<table><caption>Table 1: Geometrical and mechanical string parameters used in the simulations. R, L, $\theta$, E and $\rho$ being, respectively, the radius, length, string/tailpiece angle, Young's modulus and mass density.</caption><thead><tr><th>R (mm)</th><th>L (cm)</th><th>θ (°)</th><th>E (GPa)</th><th>ρ (kg.m<sup>-3</sup>)</th></tr></thead><tbody><tr><td>0.41</td><td>69.9</td><td>25.6</td><td>7.4</td><td>1100</td></tr></tbody></table>

The waveforms of Figs. (3a,3c,3e) show noticeable differences between the non-proportionally and proportionally damped models in terms of frequency content. In particular, there is approximately a factor of 2 between the maximum amplitudes of these models for the harp Nko. Looking at the frequency domain representations on Figs. (3b,3d,3f), two explanations are found. Firstly, amplitude differences are visible in the vicinity of harp body modes responding at the string attachment point and whose MPC is below 0.95 (indicated by colored circles placed above the curves) showing that neglecting the coupling by damping between mode combinations not verifying the separation criterion of Eq. (15) leads to an erroneous dynamic response of the system in the frequency ranges concerned. Secondly, and most importantly, when string partial coincide with complex body modes (for example mode 19 of harp Fg_2 at 754 Hz or 23 and 34 of harp Nko at 496 Hz and 733 Hz), the erroneous responses formerly addressed directly impact the amplitude of these string partials who represent the main modal contributions of the simulated signal. This leads to amplitude differences between both harp body models of 5 to 8 dB.
<figure><img src="image_4.png" /><figcaption>Figure 3: Simulated acceleration at the attachment point of the 5th string on the tailpiece for each harp of the corpus. Green full/Red dotted lines: simulations obtained with the non-proportionally/proportionally damped harp body model. Colored circles: harp body modes around which simulations show significant amplitude differences. (a) (c) (e) Waveforms; (b) (d) (f) Spectrum.</figcaption></figure>

# 4 Conclusion

The paper presents a novel modal Udwadia-Kalaba formulation based on complex modes of the dissipative subsystems. This new formulation is of particular interest in the frame of experimental substructuring, where physical and model damping matrices of substructures are usually not well estimated. Indeed, only a diagonal first order model based on eigenvectors of the state-space representation of each subsystem is required
instead of a second order model relying on a potentially fully populated modal damping matrix in case of non-proportional damping. Moreover, the assumption of proportional damping is expected to be less and less reasonable with the increase of modal density so that the modal damping matrix should be more and more populated with the increase of frequency, as long as a modal description of the substructures remains consistent.

The main limitation of the complex formulation presented here concerns substructures with defective eigenvalues. While most eigenvalues have a probability of zero for real systems, null eigenvalues associated with rigid body modes may be defective. Since the latter are frequently encountered in dynamic substructuring, their inclusion in the formulation presented in this paper would highly broaden its generality.

## Acknowledgements

This work, part of the project NGombi, was funded by the Agence Nationale de la Recherche (French National research agency), grant number ANR-19-CE27-0013-01.

## References

[1] D. de Klerk, D. J. Rixen, and S. N. Voormeeren, “General Framework for Dynamic Substructuring: History, Review and Classification of Techniques,” AIAA Journal, vol. 46, no. 5, pp. 1169–1181, May 2008, doi: 10.2514/1.33274.

[2] A. Laulusa and O. A. Bauchau, “Review of Classical Approaches for Constraint Enforcement in Multi-body Systems,” J. Comput. Nonlinear Dyn., vol. 3, no. 1, p. 011004, Nov. 2007, doi: 10.1115/1.2803257.

[3] G. A. Maggi, Principii della teoria matematica del movimento dei corpi: corso di meccanica razionale. Milano: U. Hoepli, 1896.

[4] J. I. NeMARK and N. A. Fufaev, Dynamics of nonholonomic systems, vol. 33. Providence, RI: Amer. Math. Soc., 1972.

[5] F. E. Udwadia and R. E. Kalaba, “A new perspective on constrained motion,” Proc. R. Soc. Lond. A, vol. 439, no. 1906, pp. 407–410, 1992, doi: 10.1098/rspa.1992.0158.

[6] F. E. Udwadia and R. E. Kalaba, “The Geometry of Constrained Motion,” Z. angew. Math. Mech., vol. 75, no. 8, pp. 637–640, Jan. 1995, doi: 10.1002/zamm.19950750823.

[7] F. E. Udwadia, “Equations of motion for mechanical systems: A unified approach,” Int. J. Non-Linear Mech., vol. 31, no. 6, pp. 951–958, Nov. 1996, doi: 10.1016/S0020-7462(96)00116-3.

[8] F. E. Udwadia and R. E. Kalaba, “What is the General Form of the Explicit Equations of Motion for Constrained Mechanical Systems?,” J. Appl. Mech., vol. 69, no. 3, pp. 335–339, May 2002, doi: 10.1115/1.1459071.

[9] A. Arabyan and F. Wu, “An Improved Formulation for Constrained Mechanical Systems,” Multibody Sys. Dyn., vol. 2, no. 1, pp. 49–69, Mar. 1998, doi: 10.1023/A:1009724704839.

[10] B. Jetmundsen, R. L. Bielawa, and W. G. Flannelly, “Generalized Frequency Domain Substructure Synthesis,” J. Am. Helicopter Soc., vol. 33, no. 1, pp. 55–64, Jan. 1988, doi: 10.4050/JAHS.33.1.55.

[11] T.-J. Su and J.-N. Juang, “Substructure system identification and synthesis,” J. Guid. Control Dyn., vol. 17, no. 5, pp. 1087–1095, 1994, doi: 10.2514/3.21314.

[12] M. Imregun and D. A. Robb, “Structural modification via FRF coupling using measured data,” in Proceedings of the 10th International Modal Analysis Conference, Society for Experimental Mechanics, Bethel, CT, Feb. 1992, pp. 1095–1099.
[13] J. R. Crowley, A. L. Klosterman, G. T. Rocklin, and H. Vold, “Direct structural modification using frequency response functions,” Proceedings of the 2nd International Modal Analysis Conference (IMAC II), pp. 58–65, 1984.

[14] D. de Klerk, D. J. Rixen, and J. de Jong, “The frequency based substructuring (FBS) method reformulated according to the dual domain decomposition method,” in Proceedings of the XXIV International Modal Analysis Conference, (IMAC XXIV), St Louis, MO: Society for Experimental Mechanics, Bethel, CT, 2006, pp. 1–14.

[15] R. Craig, Jr., “Coupling of substructures for dynamic analyses - An overview,” in Proceedings of the 41st Structures, Structural Dynamics, and Materials Conference, Atlanta, GA, U.S.A.: AIAA/ASME/ASCE/AHS/ASC, Apr. 2000. doi: 10.2514/6.2000-1573.

[16] R. R. Craig and M. C. C. Bampton, “Coupling of substructures for dynamic analyses.,” AIAA Journal, vol. 6, no. 7, pp. 1313–1319, 1968, doi: 10.2514/3.4741.

[17] R. H. MacNeal, “A hybrid method of component mode synthesis,” Comput. Struct., vol. 1, no. 4, pp. 581–601, Dec. 1971, doi: 10.1016/0045-7949(71)90031-9.

[18] S. Rubin, “Improved Component-Mode Representation for Structural Dynamic Analysis,” AIAA Journal, vol. 13, no. 8, pp. 995–1006, 1975, doi: 10.2514/3.60497.

[19] J. Antunes and V. Debut, “Dynamical computation of constrained flexible systems using a modal Udwadia-Kalaba formulation: Application to musical instruments,” J. Acoust. Soc. Am., vol. 141, no. 2, pp. 764–778, Feb. 2017, doi: 10.1121/1.4973534.

[20] V. Debut and J. Antunes, “Physical synthesis of six-string guitar plucks using the Udwadia-Kalaba modal formulation,” J. Acoust. Soc. Am., vol. 148, no. 2, pp. 575–587, Aug. 2020, doi: 10.1121/10.0001635.

[21] J.-T. Jiolat, C. d'Alessandro, J.-L. Le Carrou, and J. Antunes, “Toward a physical model of the clavichord,” J. Acoust. Soc. Am., vol. 150, no. 4, pp. 2350–2363, 2021.

[22] J.-T. Jiolat, J.-L. Le Carrou, and C. d'Alessandro, “Whistling in the clavichord,” J. Acoust. Soc. Am., vol. 153, pp. 338–347, Jan. 2023, doi: 10.1121/10.0016825.

[23] F. Fabre, J.-L. Le Carrou, and B. Chomette, “Physical modelling of a harp from Central Africa,” J. Acoust. Soc. Am., vol. 154, no. 4, pp. 2337–2348, Oct. 2023, doi: 10.1121/10.0021332.

[24] J. Antunes, V. Debut, L. Borsoi, X. Delaune, and P. Piteau, “A modal Udwadia-Kalaba formulation for vibro-impact modelling of continuous flexible systems with intermittent contacts,” Procedia Engineering, vol. 199, pp. 322–329, Jan. 2017, doi: 10.1016/j.proeng.2017.09.058.

[25] A. A. Shabana, “Forms of the dynamic equations,” in Computational Dynamics, John Wiley & Sons, Ltd, 2010, pp. 177–210. doi: https://doi.org/10.1002/9780470686850.ch4.

[26] D. de Falco, E. Pennestrì, and L. Vita, “Investigation of the Influence of Pseudoinverse Matrix Calculations on Multibody Dynamics Simulations by Means of the Udwadia-Kalaba Formulation,” J. Aerosp. Eng., vol. 22, no. 4, pp. 365–372, Oct. 2009, doi: 10.1061/(ASCE)0893-1321(2009)22:4(365).

[27] E. Balmés, “New results on the identification of normal modes from experimental complex modes,” Mech. Syst. Signal Process., vol. 11, no. 2, pp. 229–243, Mar. 1997, doi: 10.1006/mssp.1996.0058.

[28] T. K. Hasselman, “Modal coupling in lightly damped structures,” AIAA Journal, vol. 14, no. 11, pp. 1627–1628, 1976, doi: 10.2514/3.7259.

[29] D. Bini, B. Iannazzo, and B. Meini, Numerical solution of algebraic Riccati equations. Philadelphia: SIAM, 2011.