Dynamic substructuring

Dynamic Substructuring (DS) is an engineering tool used to model and analyse the dynamics of mechanical systems by means of its components or substructures. Using the dynamic substructuring approach one is able to analyse the dynamic behaviour of substructures separately and to later on calculate the assembled dynamics using coupling procedures. Dynamic substructuring has several advantages over the analysis of the fully assembled system: Dynamic substructuring is particularly tailored to simulation of mechanical vibrations, which has implications for many product aspects such as sound / acoustics, fatigue / durability, comfort and safety. Also, dynamic substructuring is applicable to any scale of size and frequency. It is therefore a widely used paradigm in industrial applications ranging from automotive and aerospace engineering to design of wind turbines and high-tech precision machinery.
 * Substructures can be modelled in the domain that is most appropriate, e.g. experimentally obtained substructures can be combined with numerical models.
 * Large and/or complex systems can be optimized on substructure level.
 * Numerical computation load can be reduced as solving several substructures is computationally less demanding than solving one large system.
 * Substructure models of different development groups can be shared and combined without exposing the modelling details.

History
The roots of dynamic substructuring can be found in the field of domain decomposition. In 1890 the mathematician Hermann Schwarz came up with an iterative procedure for domain decomposition which allows to solve for continuous coupled subdomains. However, many of the analytical models of coupled continuous subdomains do not have closed-form solutions, which led to discretization and approximation techniques such as the Ritz method (which is sometimes called the Rayleigh-Ritz method due to the similarity between Ritz's formulation and the Rayleigh ratio) the boundary element method (BEM) and the finite element method (FEM). These methods can be considered as "first level" domain decomposition techniques.

The finite element method proved to be the most efficient method and the invention of the microprocessor made it possible to easily solve a large variety of physical problems. In order to analyse even larger and more complex problems, methods were invented to optimize the efficiency of the discretized calculations. The first step was replacing the direct solvers by iterative solvers such as the conjugate gradient method. The lack of robustness and slow convergence of these solvers did not make them an interesting alternative in the beginning. The rise of parallel computing in the 1980s however sparked their popularity. Complex problems could now be solved by dividing the problem into subdomains, each processed by a separate processor, and solving for the interface coupling iteratively. This can be seen as a second level domain decomposition as is visualized in the figure.

The efficiency of dynamic modelling could be increased even further by reducing the complexity of the individual subdomains. This reduction of the subdomains (or substructures in the context of structural dynamics) is realized by representing substructures by means of their general responses. Expressing the separate substructures by means of their general response instead of their detailed discretization led to the so-called dynamic substructuring method. This reduction step also allowed for replacing the mathematical description of the domains by experimentally obtained information. This reduction step is also visualized by the reduction arrow in the figure.

The first dynamic substructuring methods were developed in the 1960s and were more commonly known under the name component mode synthesis (CMS). The benefits of dynamic substructuring were quickly discovered by the scientific and engineering communities and it became an important research topic in the field of structural dynamics and vibrations. Major developments followed, resulting in e.g. the classic Craig-Bampton method. The Craig-Bampton method employs static condensation (Guyan reduction) and modal truncation techniques to effectively reduce the degrees of freedom in a system.

Due to improvements in sensor and signal processing technology in the 1980s, substructuring techniques also became attractive for the experimental community. Methods dealing with structural dynamic modification were created in which coupling techniques were directly applied to measured frequency response functions (FRFs). Broad popularity of the method was gained when Jetmundsen et al. formulated the classical frequency-based substructuring (FBS) method, which laid the ground work for frequency-based dynamic substructuring. In 2006 a systematic notation was introduced by De Klerk et al. in order to simplify the difficult and elaborate notation that had been used prior. The simplification was done by means of two Boolean matrices that handle all the "bookkeeping" involved in the assembly of substructures

Domains
Dynamic substructuring can best be seen as a domain-independent toolset for assembly of component models, rather than a modelling method of its own. Generally, dynamic substructuring can be used for all domains that are well suited to simulate multiple input/multiple output behaviour. Five domains that are well suited for substructuring are:
 * Physical domain
 * Modal domain
 * Frequency domain
 * Time domain
 * State-space domain

The physical domain concerns methods that are based on (linearised) mass, damping and stiffness matrices, typically obtained from numerical FEM modelling. Popular solutions to solve the associated system of second order differential equations are the time integration schemes of Newmark and the Hilbert-Hughes-Taylor scheme. The modal domain concerns component mode synthesis (CMS) techniques such as the Craig-Bampton, Rubin and McNeal method. These methods provide efficient modal reduction bases and assembly techniques for numerical models in the physical domain. The frequency domain is more popularly known as Frequency Based Substructuring (FBS). Based on the classic formulation of Jetmundsen et al. and the reformulation of De Klerk et al., it has become the most commonly used domain for substructuring, because of the ease of expressing the differential equations of a dynamical system (by means of Frequency Response Functions, FRFs) and the convenience of implementing experimentally obtained models. The time domain refers to the recently proposed concept of Impulse Based Substructuring (IBS), which expresses the behaviour of a dynamic system using a set of Impulse Response Functions (IRFs). The state-space domain, finally, refers to methods proposed by Sjövall et al. that employ system identification techniques common to control theory.

An overview of the governing equations of the five herementioned domains is presented in the table below.

As dynamic substructuring is a domain-independent toolset, it is applicable to the dynamic equations of all domains. In order to establish substructure assembly in a particular domain, two interface conditions need to be implemented. This is explained next, followed by a few common substructuring techniques.

Interface Conditions
To establish substructuring coupling / decoupling in each of the above-mentioned domains, two conditions should be met: These are the two essential conditions that keep substructures together, hence allow to construct an assembly of multiple components. Note that the conditions are comparable with Kirchhoff's laws for electric circuits, in which case similar conditions apply to currents and voltages though/over electric components in a network; see also Mechanical–electrical analogies.
 * Coordinate compatibility, i.e. the connecting nodes of two substructures should have equal interface displacement.
 * Force equilibrium, i.e. the interface forces between connecting nodes have equal magnitude and opposing sign.

Substructure connectivity
Consider two substructures A and B as depicted in the figure. The two substructures comprise a total of six nodes; the displacements of the nodes are described by a set of Degrees of Freedom (DoFs). The DoFs of the six nodes are partitioned as follows: Note that the denotation 1, 2 and 3 indicates the function of the nodes/DoFs rather than the total amount. Let us define the sets of DoFs for the two substructures A and B in concatenated form. The displacements and applied forces are represented by the sets $$\mathbf u$$ and $$\mathbf f$$. For the purpose of substructuring, a set of interface forces $$\mathbf g$$ is introduced which only contains non-zero entries on the interface DoFs:
 * 1) DoFs of the internal nodes of substructure A;
 * 2) DoFs of the coupling nodes of substructures A and B, i.e. interface DoFs;
 * 3) DoFs of the internal nodes of substructure B.

\mathbf u \triangleq \begin{bmatrix} \mathbf u_1^A \\ \mathbf u_2^A \\ \mathbf u_2^B \\ \mathbf u_3^B \\ \end{bmatrix} \quad \mathbf f \triangleq \begin{bmatrix} \mathbf f_1^A \\ \mathbf f_2^A \\ \mathbf f_2^B \\ \mathbf f_3^B \\ \end{bmatrix} \quad \mathbf g \triangleq \begin{bmatrix} \mathbf 0 \\ \mathbf g_2^A \\ \mathbf g_2^B \\ \mathbf 0 \\ \end{bmatrix} \quad \mathbf u, \mathbf f,\mathbf g \in \R^n $$ The relation between dynamic displacements $$\mathbf u$$ and applied forces $$\mathbf f$$ of the uncoupled problem is governed by a particular dynamic equation, such as presented in the table above. The uncoupled equations of motion are augmented by extra terms/equations for compatibility and equilibrium, as discussed next.

Compatibility
The compatibility condition requires that the interface DoFs have the same sign and value at both sides of the interface: $$ \mathbf u_2^A = \mathbf u_2^B $$. This condition can be expressed using a so-called signed Boolean matrix, denoted by $$\mathbf B $$. For the given example this can be expressed as:

\mathbf B \mathbf u = \mathbf 0 \quad \Rightarrow \quad \mathbf u_2^B - \mathbf u_2^A = \mathbf 0 \quad \Rightarrow \quad \mathbf B \triangleq \begin{bmatrix} \mathbf 0 & -\mathbf I & \mathbf I & \mathbf 0 \end{bmatrix} $$

In some cases the interface nodes of the substructures are non-conforming, e.g. when two substructures are meshed separately. In such cases a non-Boolean matrix $$ \mathbf B $$ has to be used in order to enforce a weak interface compatibility.

A second form in which the compatibility condition can be expressed is by means of coordinate substitution by a set of generalised coordinates $$ \mathbf q $$. The set $$ \mathbf q $$ contains the unique coordinates that remain after assembly of the substructures. Every matching pair of interface DoFs is described by a single generalised coordinate, which means that the compatibility condition is automatically enforced. Expressing $$ \mathbf u $$ using $$ \mathbf q $$ gives:



\mathbf u = \mathbf L \mathbf q \quad \Rightarrow \quad \begin{cases} \mathbf u_1^A = \mathbf q_1 \\ \mathbf u_2^A = \mathbf q_2 \\ \mathbf u_2^B = \mathbf q_2 \\ \mathbf u_3^B = \mathbf q_3 \\ \end{cases} \quad \Rightarrow \quad \mathbf L \triangleq \begin{bmatrix} \mathbf I & \mathbf 0 & \mathbf 0 \\ \mathbf 0 & \mathbf I & \mathbf 0 \\ \mathbf 0 & \mathbf I & \mathbf 0 \\ \mathbf 0 & \mathbf 0 & \mathbf I \\ \end{bmatrix} $$

Matrix $$ \mathbf L $$ is referred to as the Boolean localisation matrix. A useful relation between matrix $$ \mathbf B $$ and $$ \mathbf L $$ can be exposed by noting that compatibility should hold for any set of physical coordinates $$\mathbf u$$ expressed by $$\mathbf q$$. Indeed, substituting $$\mathbf u = \mathbf{Lq}$$ in the equation $$\mathbf{Bu} = \mathbf 0$$:

\mathbf{Bu} = \mathbf B \mathbf L \mathbf q = \mathbf 0 \quad \forall \mathbf q $$

Hence $$ \mathbf L $$ represents the nullspace of $$ \mathbf B $$:

\begin{cases} \mathbf L = \text{null}(\mathbf B)\\ \mathbf B^T = \text{null}(\mathbf L^T) \end{cases} $$

This means in practice that one only needs to define $$\mathbf B$$ or $$\mathbf L$$; the other Boolean matrix is calculated using the nullspace property.

Equilibrium
The second condition that has to be satisfied for substructure assembly is the force equilibrium for matching interface forces $$\mathbf g$$. For the current example, this condition can be written as $$ \mathbf g_2^A = -\mathbf g_2^B $$. Similar to the compatibility equation, the force equilibrium condition can be expressed using a Boolean matrix. Use is made of the transpose of the Boolean localisation matrix $$ \mathbf L $$ that was introduced to write compatibility:

\mathbf L^T \mathbf g = \mathbf 0 \quad \Rightarrow \quad \begin{cases} \mathbf g_1^A = \mathbf 0 \\ \mathbf g_2^A + \mathbf g_2^B = \mathbf 0 \\ \mathbf g_3^B = \mathbf 0 \end{cases} $$

The equations for $$\mathbf g_1$$ and $$\mathbf g_3$$ state that the interface forces on internal nodes are zero, hence not present. The equation for $$\mathbf g_2$$ correctly establishes the force equilibrium between a matching pair of interface DoFs according to Newton's third law.

A second notation in which the equilibrium condition can be expressed is by introducing a set of Lagrange multipliers $$\boldsymbol \lambda $$. The substitution of these Lagrange multipliers is possible as $$ \mathbf g_2^A $$ and $$ \mathbf g_2^B $$ differ only in sign, not in value. Using again the signed Boolean matrix $$ \mathbf B $$:



\mathbf g = -\mathbf B^T \boldsymbol{\lambda} \quad \Rightarrow \quad \begin{cases} \mathbf g_1^A = \mathbf 0 \\ \mathbf g_2^A = \boldsymbol{\lambda}\\ \mathbf g_2^B = -\boldsymbol{\lambda}\\ \mathbf g_3^B =\mathbf 0 \\ \end{cases} $$

The set $$\boldsymbol{\lambda}$$ defines the intensity of the interface forces $$\mathbf g_2$$. Each Lagrange multiplier represents the magnitude of two matching interface forces in the assembly. By defining the interface forces $$\mathbf g$$ using Lagrange multipliers $$\boldsymbol \lambda $$, force equilibrium is automatically satisfied. This can be seen by substituting $$\mathbf g = -\mathbf{B^T} \boldsymbol{\lambda}$$ into the first equilibrium equation:

\mathbf L^T \mathbf g = -\mathbf L^T \mathbf B^T \boldsymbol{\lambda} = \mathbf 0 \quad \forall \mathbf \mathbf g $$

Again, the nullspace property of the Boolean matrices is used here, namely: $$ \mathbf{L^T}\mathbf{B^T} = \mathbf 0 $$.

The two conditions as presented above can be applied to establish coupling / decoupling in a myriad of domains and are thus independent of variables such as time, frequency, mode, etc. Some implementations of the interface conditions for the most common domains of substructuring are presented below.

Substructuring in the physical domain
The physical domain is the domain that has the most straightforward physical interpretation. For each discrete linearised dynamic system one is able to write an equilibrium between the externally applied forces and the internal forces originating from intrinsic inertia, viscous damping and elasticity. This relation is governed by one of the most elementary formulas in structural vibrations:


 * $$\mathbf{M} \mathbf \ddot{u}(t)+\mathbf C \mathbf \dot{u}(t)+ \mathbf K \mathbf u(t) = \mathbf f(t)$$

$$\mathbf M, \mathbf C, \mathbf K$$ represent the mass, damping and stiffness matrix of the system. These matrices are often obtained from finite element modelling (FEM), and are referred to as the numerical model of the structure. Furthermore, $$\mathbf u$$ represents the DoFs and $$\mathbf f$$ the force vector which are dependent on time $$(t)$$. This dependency is omitted in the following equations in order to improve readability.

Coupling in the physical domain
Coupling of $$n$$ substructures in the physical domain first requires writing the uncoupled equations of motion of the $$n$$ substructures in block diagonal form:


 * $$\mathbf M \triangleq \text{diag}(\mathbf M^{(1)},

\dots, \mathbf M^{(n)}) = \begin{bmatrix} \mathbf M^{(1)} & . & . \\ . & \ddots & . \\ . & . & \mathbf M^{(n)} \\ \end{bmatrix}$$
 * $$\mathbf C \triangleq \text{diag}(\mathbf C^{(1)},\dots,\mathbf C^{(n)}) \quad \mathbf K \triangleq \text{diag}(\mathbf K^{(1)},\dots,\mathbf K^{(n)})$$
 * $$\mathbf u \triangleq

\begin{bmatrix} \mathbf u^{(1)} \\ \vdots \\ \mathbf u^{(n)} \end{bmatrix}, \quad \mathbf f \triangleq \begin{bmatrix} \mathbf f^{(1)} \\ \vdots \\ \mathbf f^{(n)} \end{bmatrix}, \quad \mathbf g \triangleq \begin{bmatrix} \mathbf g^{(1)}\\ \vdots \\ \mathbf g^{(n)} \end{bmatrix}$$

Next, two assembly approaches can be distinguished: primal and dual assembly.

Primal assembly
For primal assembly, a unique set of degrees of freedom $$\mathbf q$$ is defined in order to satisfy compatibility, $$ \mathbf u = \mathbf {Lq} $$. Furthermore, a second equation is added to enforce interface force equilibrium. This results in the following coupled dynamic equilibrium equations:


 * $$\begin{cases}

\mathbf{ML\ddot{q} + CL\dot{q}+KLq = f + g}\\ \mathbf{L^Tg = 0} \end{cases}$$

Pre-multiplying the first equation by $$\mathbf L^T$$ and noting that $$\mathbf{L^Tg}= \mathbf{0}$$, the primal assembly reduces to:


 * $$\mathbf{\tilde{M}\ddot{q}+\tilde{C}\dot{q}+\tilde{K}q = L^T f} \quad

\begin{cases} \mathbf{\tilde{M} = L^TML} \\ \mathbf{\tilde{C} = L^TCL} \\ \mathbf{\tilde{K} = L^TKL} \end{cases}$$

The primally assembled system matrices can be used for a transient simulation by any standard time stepping algorithm. Note that the primal assembly technique is analogue to assembly of super-elements in finite element methods.

Dual assembly
In the dual assembly formulation the global set of DoFs is retained and an assembly is made by a priori satisfying the equilibrium condition $$\mathbf{g = -B^T \boldsymbol{\lambda}}$$. Again, the Lagrange multipliers represent the interface forces connecting the DoFs at the interface. As these are unknowns, they are moved to the left-hand side of the equation. In order to satisfy compatibility, a second equation is added to the system, now operating on the displacements:


 * $$\begin{cases}

\mathbf{M\ddot{u} + C \dot{u} + K u + B^T \boldsymbol{\lambda} = f}\\ \mathbf{Bu = 0} \end{cases}$$

The dually assembled system can be written in matrix form as:


 * $$\mathbf{\begin{bmatrix}

\mathbf M & \mathbf 0 \\ \mathbf  0 &\mathbf  0 \end{bmatrix}} \begin{bmatrix} \mathbf \ddot{u} \\ \boldsymbol{\lambda} \end{bmatrix} + \begin{bmatrix} \mathbf C & \mathbf 0 \\ \mathbf  0 & \mathbf  0 \end{bmatrix} \begin{bmatrix} \mathbf \dot{u} \\ \boldsymbol{\lambda} \end{bmatrix} + \begin{bmatrix} \mathbf K & \mathbf B^T \\ \mathbf B & \mathbf 0 \end{bmatrix} \begin{bmatrix} \mathbf u \\ \boldsymbol{\lambda} \end{bmatrix} = \begin{bmatrix} \mathbf f \\ \mathbf 0 \end{bmatrix}$$

This dually assembled system can also be used in a transient simulation by means of a standard time stepping algorithm.

Substructuring in the frequency domain
In order to write out the equations for frequency based substructuring (FBS), the dynamic equilibrium first has to be put in the frequency domain. Starting with the dynamic equilibrium in the physical domain:
 * $$\mathbf{M} \mathbf \ddot{u}(t)+\mathbf C \mathbf \dot{u}(t)+ \mathbf K \mathbf u(t) = \mathbf f(t)$$

Taking the Fourier transform of this equation gives the dynamic equilibrium in the frequency domain:



\mathbf Z(\omega) \mathbf u(\omega) = \mathbf f(\omega) \quad \text{with} \quad \mathbf Z(\omega) = \bigr[-\omega^2 \mathbf M+j\omega \mathbf C+\mathbf K\bigr] $$

Matrix $$\mathbf Z(\omega)$$ is referred to as the dynamic stiffness matrix. This matrix consists of the complex-valued frequency-dependent functions that describe the force required to generate a unit harmonic displacement at a certain DoF. The inverse of the matrix $$\mathbf Z(\omega)$$ is defined as $$\mathbf Y(\omega)\triangleq(\mathbf Z(\omega))^{-1}$$ and yields the more intuitive admittance notation:



\mathbf u(\omega) = \mathbf Y(\omega) \mathbf f(\omega) $$

The receptance matrix $$\mathbf Y(\omega)$$ contains the frequency response functions (FRFs) of the structure which describe the displacement response to a unit input force. Other variants of the receptance matrix are the mobility and accelerance matrix, which respectively describe the velocity and acceleration response. The elements of the dynamic stiffness (or impedance in general) and receptance (or admittance in general) matrix are defined as follows:



\begin{align} \text{Impedance:} \quad Z_{ij}(\omega)= \frac{f_i(\omega)}{u_j(\omega)} \quad u_{k\neq j}=0\\ \text{Admittance:} \quad Y_{ij}(\omega)= \frac{u_i(\omega)}{f_j(\omega)}\quad f_{k\neq j}=0 \end{align} $$

Coupling in the frequency domain
In order to couple two substructures in the frequency domain, use is made of the admittance and impedance matrices of both substructures. Using the definition of substructures A and B as introduced previously, the following impedance and admittance matrices are defined (note that the frequency dependency $$(\omega) $$ is omitted from the terms to improve readability):


 * $$\begin{array}{ll}

\mathbf Z^A \triangleq \begin{bmatrix} \mathbf Z_{11}^A & \mathbf Z_{12}^A \\ \mathbf Z_{21}^A & \mathbf Z_{22}^A \\ \end{bmatrix}

&

\mathbf Z^B \triangleq \begin{bmatrix} \mathbf Z_{22}^B & \mathbf Z_{23}^B \\ \mathbf Z_{32}^B & \mathbf Z_{33}^B \\ \end{bmatrix}

\\[18pt]

\mathbf Y^A \triangleq \begin{bmatrix} \mathbf Y_{11}^A & \mathbf Y_{12}^A \\ \mathbf Y_{21}^A & \mathbf Y_{22}^A \\ \end{bmatrix}

& \mathbf Y^B \triangleq \begin{bmatrix} \mathbf Y_{22}^B & \mathbf Y_{23}^B \\ \mathbf Y_{32}^B & \mathbf Y_{33}^B \\ \end{bmatrix}

\end{array} $$

The two admittance and impedance matrices can be put in block diagonal form in order to align with the global set of DoFs $$\mathbf u $$:


 * $$\begin{align}

\mathbf Z \triangleq &\begin{bmatrix} \mathbf Z^A &\mathbf 0 \\ \mathbf0 & \mathbf Z^B \end{bmatrix} \triangleq \begin{bmatrix} \mathbf Z_{11}^A & \mathbf Z_{12}^A & \mathbf 0 & \mathbf 0 \\ \mathbf Z_{21}^A & \mathbf Z_{22}^A & \mathbf 0 & \mathbf 0 \\ \mathbf 0 & \mathbf 0 & \mathbf Z_{22}^B & \mathbf Z_{23}^B \\ \mathbf 0 & \mathbf0 & \mathbf Z_{32}^B & \mathbf Z_{33}^B \\ \end{bmatrix} \\[6pt]

\mathbf Y \triangleq &\begin{bmatrix} \mathbf Y^A & \mathbf0 \\ \mathbf 0 & \mathbf Y^B \end{bmatrix} \triangleq \begin{bmatrix} \mathbf Y_{11}^A & \mathbf Y_{12}^A & \mathbf 0 & \mathbf 0 \\ \mathbf Y_{21}^A & \mathbf Y_{22}^A &\mathbf 0 & \mathbf 0 \\ \mathbf 0 & \mathbf 0 & \mathbf Y_{22}^B & \mathbf Y_{23}^B \\ \mathbf 0 & \mathbf 0 & \mathbf Y_{32}^B & \mathbf Y_{33}^B \\ \end{bmatrix} \end{align}$$

The off-diagonal zero terms show that at this point no coupling is present between the two substructures. In order to create this coupling, use can be made of the primal or dual assembly method. Both assembly methods make use of the dynamic equations as was defined before:


 * $$\begin{align}

\mathbf Z \mathbf u & = \mathbf f + \mathbf g\\ \mathbf u & = \mathbf Y(\mathbf f+ \mathbf g) \\ \end{align}$$

In these equations $$\mathbf g$$ is again used to define the set of interface forces, which are yet unknown.

Primal assembly
In order to obtain the primal system of equations, a unique set of coordinates $$ \mathbf q $$ is defined: $$\mathbf u = \mathbf{Lq}$$. By definition of an appropriate Boolean localisation matrix $$\mathbf L$$, a unique set of DoFs remains for which the compatibility condition is satisfied a priori (compatibility condition). In order to satisfy the equilibrium condition a second equation is added to the equations of motion:


 * $$\begin{cases}

\mathbf Z \mathbf L \mathbf q = \mathbf f + \mathbf g \\ \mathbf L^T \mathbf g = \mathbf 0 \end{cases}$$

Pre-multiplying the first equation with $$\mathbf L^T$$ yields the notation of the assembled equations of motion for the generalised coordinates $$\mathbf q$$:


 * $$\mathbf \tilde{Z}\mathbf q = \mathbf \tilde{f}

\quad \text{with} \quad \begin{cases} \mathbf \tilde{Z} = \mathbf L^T\mathbf Z \mathbf L \\ \mathbf \tilde{f} = \mathbf L^T \mathbf f \end{cases}$$

This result can be rewritten in admittance form as:


 * $$\mathbf q = \mathbf \tilde{Y}\mathbf \tilde{f}

\quad \text{with} \quad \mathbf \tilde{Y}\ = (\mathbf \tilde{Z})^{-1}$$

This last result gives access to the generalised responses as a result of the generalised applied forces $$\mathbf \tilde{f}$$, namely by inverting the primally assembled impedance matrix.

The primal assembly procedure is mainly of interest when one has access to the dynamics in impedance form, e.g. from finite element modelling. When one only has access to the dynamics in admittance notation, the dual formulation is a more suitable approach.

Dual assembly
A dually assembled system starts with the system written in the admittance notation. For a dually assembled system the force equilibrium condition is satisfied a priori by substituting Lagrange multipliers $$\boldsymbol{\lambda}$$ for the interface forces: $$\mathbf g = - \mathbf B^T\boldsymbol{\lambda}$$. The compatibility condition is enforced by adding an additional equation:

\begin{cases} \mathbf u = \mathbf Y(\mathbf f - \mathbf B^T \boldsymbol{\lambda})\\ \mathbf B \mathbf u = \mathbf 0 \end{cases} $$

Substituting the first line in the second and solving for $$ \boldsymbol{\lambda} $$ gives:

\boldsymbol{\lambda} = (\mathbf B \mathbf Y \mathbf B^T)^{-1}\mathbf B \mathbf Y \mathbf f $$

The term $$ \mathbf B \mathbf Y \mathbf f $$ represents the incompatibility caused by the uncoupled responses of the substructures to the applied forces $$\mathbf f $$. By multiplying the incompatibility with the combined interface stiffness, i.e. $$ (\mathbf B \mathbf  Y \mathbf B^T)^{-1}$$, the forces $$ \boldsymbol{\lambda} $$ that keep the substructures together are determined. The coupled response is obtained by substituting the calculated $$ \boldsymbol{\lambda} $$ back into the original equation:

\begin{align} &\mathbf u = \mathbf Y \mathbf f - \mathbf Y \mathbf B^T(\mathbf{BYB^T})^{-1}\mathbf{BYf} \\[3pt] &\mathbf u = \mathbf \tilde{Y} \mathbf f \quad \text{with:} \quad \mathbf \tilde{Y} = \big(\mathbf I - \mathbf Y \mathbf B^T(\mathbf{BYB^T})^{-1}\mathbf{B}\big)\mathbf Y \end{align} $$

This coupling method is referred to as the Lagrange-multiplier frequency-based substructuring (LM-FBS) method. The LM-FBS method allows for quick and easy assembling of an arbitrary number of substructures in a systematic fashion. Note that the result $$ \mathbf{\tilde{Y}} $$ is theoretically the same as was obtained above by application of primal assembly.



Decoupling in the frequency domain
In addition to coupling of substructures, one is also able to decouple substructures from assemblies. Using the plus sign as a substructure coupling operator, the coupling procedure could simply be described as AB = A + B. Using a similar notation, decoupling could be formulated as AB - B = A. Decoupling procedures are often required to remove substructures that were added for measurement purposes, e.g. to fix the structure. Similar to coupling, a primal and dual formulation exists for decoupling procedures.

Primal disassembly
As a result of the primal coupling, the impedance matrix of the assembled system can be written as follows:

\mathbf Z^{AB}= \begin{bmatrix} \mathbf Z_{11}^{AB} & \mathbf Z_{12}^{AB} & \mathbf Z_{13}^{AB} \\ \mathbf Z_{21}^{AB} & \mathbf Z_{22}^{AB} & \mathbf Z_{23}^{AB} \\ \mathbf Z_{31}^{AB} & \mathbf Z_{32}^{AB} & \mathbf Z_{33}^{AB} \end{bmatrix} = \begin{bmatrix} \mathbf Z_{11}^A & \mathbf Z_{12}^A & \mathbf 0 \\ \mathbf Z_{21}^A & \mathbf Z_{22}^A+ \mathbf Z_{22}^B & \mathbf Z_{23}^B\\ \mathbf 0 & \mathbf Z_{32}^B & \mathbf Z_{33}^B \end{bmatrix} $$

Using this relation, the following trivial subtraction operation would suffice for the decoupling of the substructure B from assembly AB:



\begin{bmatrix} \mathbf Z_{11}^A & \mathbf Z_{12}^A & \mathbf 0 \\ \mathbf Z_{21}^A & \mathbf Z_{22}^A+\mathbf Z_{22}^B & \mathbf Z_{23}^B\\ \mathbf 0 & \mathbf Z_{32}^B & \mathbf Z_{33}^B \end{bmatrix}- \begin{bmatrix} \mathbf 0 &\mathbf 0 &\mathbf  0 \\ \mathbf 0 & \mathbf Z_{22}^B & \mathbf Z_{23}^B \\ \mathbf 0 & \mathbf Z_{32}^B & \mathbf Z_{33}^B \end{bmatrix} = \begin{bmatrix} \mathbf Z_{11}^A & \mathbf Z_{12}^A & \mathbf 0 \\ \mathbf Z_{21}^A & \mathbf Z_{22}^A & \mathbf0 \\ \mathbf0 & \mathbf0 & \mathbf0 \end{bmatrix} $$ By placing the impedance of AB and B in block-diagonal form, with a minus sign for the impedance of B to account for the subtraction operation, the same equation that was used for primal coupling can now be used to perform the primal decoupling procedures.


 * $$\mathbf Z \triangleq

\begin{bmatrix} \mathbf Z^{AB} & \mathbf0 \\ \mathbf0 & -\mathbf Z^{B} \end{bmatrix} \quad \Rightarrow \quad \mathbf \tilde{Z} = \mathbf L^T\mathbf Z \mathbf L $$

with:


 * $$\mathbf L = \begin{bmatrix}

\mathbf I & \mathbf0 & \mathbf0 \\ \mathbf0 & \mathbf I & \mathbf0 \\ \mathbf0 & \mathbf0 & \mathbf I \\ \mathbf0 & \mathbf I & \mathbf0 \\ \mathbf0 & \mathbf0 & \mathbf I \end{bmatrix}$$

The primal disassembly can thus be understood as the assembly of structure AB with the negative impedance of substructure B. A limitation of the primal disassembly is that all DoF of the substructure that is to be decoupled have to be exactly represented in the assembled situation. For numerical decoupling situations this should not pose any problems, however for experimental cases this can be troublesome. A solution to this problem can be found in the dual disassembly.

Dual disassembly
Similar to the dual assembly, the dual disassembly approaches the decoupling problem using the admittance matrices. Decoupling in the dual domain means finding a force that ensures compatibility, yet acts in the opposite direction. This newly found force would then counteract the force that is applied to the assembly due to the dynamics of substructure B. Writing this out in equations of motion:



\begin{cases} \mathbf u^{AB} = \mathbf Y^{AB} \mathbf f + \mathbf Y^{AB}\mathbf g \\ \mathbf u^B = -\mathbf Y^B \mathbf g \end{cases} $$

In order to write the dynamics of both systems in one equation, using the LM-FBS assembly notation, the following matrices are defined:


 * $$\begin{align}

\mathbf Y \triangleq &\begin{bmatrix} \mathbf Y^{AB} & \mathbf0 \\ \mathbf0 & -\mathbf Y^{B} \end{bmatrix} \\[3pt] \mathbf u = &\begin{bmatrix} \mathbf u_1^{AB} \\ \mathbf u_2^{AB} \\ \mathbf u_3^{AB} \\ \mathbf u_2^{B} \\ \mathbf u_3^{B} \end{bmatrix}; \quad \mathbf f \triangleq \begin{bmatrix} \mathbf f_1^{AB}\\ \mathbf f_2^{AB}\\ \mathbf0 \\ \mathbf0 \\ \mathbf0 \end{bmatrix}; \quad \mathbf g \triangleq \begin{bmatrix} \mathbf0 \\ \mathbf g_2^{AB}\\ \mathbf0 \\ \mathbf g_2^B \\ \mathbf0 \end{bmatrix}

\end{align} $$

In order to enforce compatibility, a similar approach is used as for the assembly task. Defining a $$\mathbf B$$-matrix to enforce compatibility:


 * $$\mathbf B = \begin{bmatrix}

\mathbf0 & -\mathbf I & \mathbf0 & \mathbf I & \mathbf0 \end{bmatrix} $$

Using this notation, the disassembly procedure can be performed using exactly the same equation as was used for the dual assembly:



\begin{cases} \mathbf u = \mathbf Y(\mathbf f - \mathbf B^T \boldsymbol{\lambda})\\ \mathbf B \mathbf u = \mathbf0 \end{cases} $$

This means that coupling and decoupling procedures using LM-FBS require identical steps, the only difference being the manner in which the global admittance matrix is defined. Indeed, the substructures to couple appear with a plus sign, whereas decoupled structures carry a minus sign:


 * $$\begin{align}

\text{coupling} \quad \mathbf Y & \triangleq \begin{bmatrix} \mathbf Y_A & \mathbf 0 \\ \mathbf 0 & \mathbf Y_B \end{bmatrix} \\ \text{decoupling} \quad \mathbf Y & \triangleq \begin{bmatrix} \mathbf Y_{AB} & \mathbf 0 \\ \mathbf 0 & -\mathbf Y_B \end{bmatrix} \end{align} $$ More advanced decoupling techniques use the fact that internal points of substructure B appear in both the admittances of AB and B, hence can be used to enhance the decoupling process. Such techniques are described in.