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.
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[1] (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.[2] 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.[3] 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.[4]
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[5]
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.[5] Five domains that are well suited for substructuring are:
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[6] and the Hilbert-Hughes-Taylor scheme.[7] 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.[8] 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),[9] 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.[10] 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.
Physical domain | f(t)=M\ddot{u | M,C,K | ||||
Modal Domain |
=
| Mm,Cm,Km \boldsymbolη | ||||
Frequency Domain | \begin{align} f(\omega)=Z(\omega)u(\omega)\\ u(\omega)=Y(\omega)f(\omega) \end{align} | Z(\omega) Y(\omega) | ||||
Time domain | u(t)=Y(t)*f(t) | Y(t) | ||||
State-space domain |
=Ax(t)+Bv(t)\\ y(t)=Cx(t)+Dv(t)\\ \end{cases} | A,B,C,D x v y |
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.
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
u
f
g
u\triangleq
A | |
\begin{bmatrix} u | |
1 |
A | |
\\ u | |
2 |
B | |
\\ u | |
2 |
B | |
\\ u | |
3 |
\\ \end{bmatrix} f\triangleq
A | |
\begin{bmatrix} f | |
1 |
A | |
\\ f | |
2 |
B | |
\\ f | |
2 |
B | |
\\ f | |
3 |
\\ \end{bmatrix} g\triangleq
A | |
\begin{bmatrix} 0\\ g | |
2 |
B | |
\\ g | |
2 |
\\ 0\\ \end{bmatrix} u,f,g\in\Rn
u
f
The compatibility condition requires that the interface DoFs have the same sign and value at both sides of the interface:
A | |
u | |
2 |
=
B | |
u | |
2 |
B
Bu=0 ⇒
B | |
u | |
2 |
-
A | |
u | |
2 |
=0 ⇒ B\triangleq\begin{bmatrix} 0&-I&I&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
B
A second form in which the compatibility condition can be expressed is by means of coordinate substitution by a set of generalised coordinates
q
q
u
q
u=Lq ⇒
A | |
\begin{cases} u | |
1 |
=q1
A | |
\\ u | |
2 |
=q2
B | |
\\ u | |
2 |
=q2
B | |
\\ u | |
3 |
=q3\\ \end{cases} ⇒ L\triangleq\begin{bmatrix} I&0&0\\ 0&I&0\\ 0&I&0\\ 0&0&I\\ \end{bmatrix}
Matrix
L
B
L
u
q
u=Lq
Bu=0
Bu=BLq=0 \forallq
Hence
L
B
\begin{cases} L=null(B)\\ BT=null(LT) \end{cases}
This means in practice that one only needs to define
B
L
The second condition that has to be satisfied for substructure assembly is the force equilibrium for matching interface forces
g
A | |
g | |
2 |
=
B | |
-g | |
2 |
L
LTg=0 ⇒
A | |
\begin{cases} g | |
1 |
=
A | |
0\\ g | |
2 |
+
B | |
g | |
2 |
=
B | |
0\\ g | |
3 |
=0 \end{cases}
The equations for
g1
g3
g2
A second notation in which the equilibrium condition can be expressed is by introducing a set of Lagrange multipliers
\boldsymbolλ
A | |
g | |
2 |
B | |
g | |
2 |
B
g=-BT\boldsymbol{λ} ⇒
A | |
\begin{cases} g | |
1 |
=
A | |
0\\ g | |
2 |
=
B | |
\boldsymbol{λ}\\ g | |
2 |
=
B | |
-\boldsymbol{λ}\\ g | |
3 |
=0\\ \end{cases}
The set
\boldsymbol{λ}
g2
g
\boldsymbolλ
g=
-BT |
\boldsymbol{λ}
LTg=-LTBT\boldsymbol{λ}=0 \forallg
Again, the nullspace property of the Boolean matrices is used here, namely:
LTBT |
=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.
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:
M
|
f(t)
M,C,K
u
f
(t)
Coupling of
n
n
M\triangleqdiag(M(1), ...,M(n)) = \begin{bmatrix} M(1)&.&.\\ .&\ddots&.\\ .&.&M(n)\\ \end{bmatrix}
C\triangleqdiag(C(1),...,C(n)) K\triangleqdiag(K(1),...,K(n))
u\triangleq\begin{bmatrix} u(1)\ \vdots\ u(n)\end{bmatrix}, f\triangleq\begin{bmatrix} f(1)\ \vdots\ f(n)\end{bmatrix}, g\triangleq\begin{bmatrix} g(1)\ \vdots\ g(n)\end{bmatrix}
Next, two assembly approaches can be distinguished: primal and dual assembly.
For primal assembly, a unique set of degrees of freedom
q
u=Lq
\begin{cases} ML\ddot{q+CL
q |
+KLq=f+
g}\\ LTg=0 \end{cases} |
Pre-multiplying the first equation by
LT
LTg= |
0
\tilde{M\ddot{q}+\tilde{C}q |
+\tilde{K}q=LTf} \begin{cases} \tilde{M=LTML}\\ \tilde{C=LTCL}\\ \tilde{K=LTKL} \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.
In the dual assembly formulation the global set of DoFs is retained and an assembly is made by a priori satisfying the equilibrium condition
g=-BT\boldsymbol{λ |
\begin{cases} M\ddot{u+C
u |
+Ku+BT\boldsymbol{λ}=f}\\ Bu=0 \end{cases}
The dually assembled system can be written in matrix form as:
\begin{bmatrix M&0\ 0&0 \end{bmatrix}} \begin{bmatrix} \ddot{u}\ \boldsymbol{λ} \end{bmatrix} + \begin{bmatrix} C&0\ 0&
|
BT\ B&0 \end{bmatrix} \begin{bmatrix} u\ \boldsymbol{λ} \end{bmatrix} = \begin{bmatrix} f\ 0 \end{bmatrix}
This dually assembled system can also be used in a transient simulation by means of a standard time stepping algorithm.
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:
M
|
f(t)
Taking the Fourier transform of this equation gives the dynamic equilibrium in the frequency domain:
Z(\omega)u(\omega)=f(\omega) with Z(\omega)=r[-\omega2M+j\omegaC+Kr]
Matrix
Z(\omega)
Z(\omega)
Y(\omega)\triangleq(Z(\omega))-1
u(\omega)=Y(\omega)f(\omega)
The receptance matrix
Y(\omega)
\begin{align} Impedance: Zij(\omega)=
fi(\omega) | |
uj(\omega) |
uk ≠ =0\\ Admittance: Yij(\omega)=
ui(\omega) | |
fj(\omega) |
fk ≠ =0 \end{align}
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)
\begin{array}{ll} ZA\triangleq
A | |
\begin{bmatrix} Z | |
11 |
&
A | |
Z | |
12 |
A | |
\\ Z | |
21 |
&
A | |
Z | |
22 |
\\ \end{bmatrix} & ZB\triangleq
B | |
\begin{bmatrix} Z | |
22 |
&
B | |
Z | |
23 |
B | |
\\ Z | |
32 |
&
B | |
Z | |
33 |
\\ \end{bmatrix} \\[18pt] YA\triangleq
A | |
\begin{bmatrix} Y | |
11 |
&
A | |
Y | |
12 |
A | |
\\ Y | |
21 |
&
A | |
Y | |
22 |
\\ \end{bmatrix} & YB\triangleq
B | |
\begin{bmatrix} Y | |
22 |
&
B | |
Y | |
23 |
B | |
\\ Y | |
32 |
&
B | |
Y | |
33 |
\\ \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
u
\begin{align} Z\triangleq &\begin{bmatrix} ZA&0\\ 0&
A | |
Z | |
11 |
&
A | |
Z | |
12 |
&0&
A | |
0\\ Z | |
21 |
&
A | |
Z | |
22 |
&0&0\\ 0&0&
B | |
Z | |
22 |
&
B | |
Z | |
23 |
\\ 0&0&
B | |
Z | |
32 |
&
B | |
Z | |
33 |
\\ \end{bmatrix}\\[6pt] Y\triangleq &\begin{bmatrix} YA&0\\ 0&
A | |
Y | |
11 |
&
A | |
Y | |
12 |
&0&
A | |
0\\ Y | |
21 |
&
A | |
Y | |
22 |
&0&0\\ 0&0&
B | |
Y | |
22 |
&
B | |
Y | |
23 |
\\ 0&0&
B | |
Y | |
32 |
&
B | |
Y | |
33 |
\\ \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} Zu&=f+g\\ u&=Y(f+g)\\ \end{align}
In these equations
g
In order to obtain the primal system of equations, a unique set of coordinates
q
u=Lq
L
\begin{cases} ZLq=f+g\\ LTg=0 \end{cases}
Pre-multiplying the first equation with
LT
q
\tilde{Z}q=\tilde{f} with \begin{cases} \tilde{Z}=LTZL\\ \tilde{f}=LTf \end{cases}
This result can be rewritten in admittance form as:
q=\tilde{Y}\tilde{f} with \tilde{Y} =(\tilde{Z})-1
This last result gives access to the generalised responses as a result of the generalised applied forces
\tilde{f}
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,[14] the dual formulation is a more suitable approach.
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{λ}
g=-BT\boldsymbol{λ}
\begin{cases} u=Y(f-BT\boldsymbol{λ})\\ Bu=0 \end{cases}
Substituting the first line in the second and solving for
\boldsymbol{λ}
\boldsymbol{λ}=(BYBT)-1BYf
The term
BYf
f
(BYBT)-1
\boldsymbol{λ}
\boldsymbol{λ}
\begin{align} &u=Yf-
| ||
YB |
-1BYf\\[3pt] &u=\tilde{Y}f with: \tilde{Y}=(I-
| ||
YB |
-1B)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
\tilde{Y
In addition to coupling of substructures, one is also able to decouple substructures from assemblies.[15] [16] [17] [18] 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.
As a result of the primal coupling, the impedance matrix of the assembled system can be written as follows:
ZAB
AB | |
= \begin{bmatrix} Z | |
11 |
&
AB | |
Z | |
12 |
&
AB | |
Z | |
13 |
AB | |
\\ Z | |
21 |
&
AB | |
Z | |
22 |
&
AB | |
Z | |
23 |
AB | |
\\ Z | |
31 |
&
AB | |
Z | |
32 |
&
AB | |
Z | |
33 |
\end{bmatrix} =
A | |
\begin{bmatrix} Z | |
11 |
&
A | |
Z | |
12 |
&
A | |
0\\ Z | |
21 |
&
A+ | |
Z | |
22 |
B | |
Z | |
22 |
&
B\\ 0& | |
Z | |
23 |
B | |
Z | |
32 |
&
B \end{bmatrix} | |
Z | |
33 |
Using this relation, the following trivial subtraction operation would suffice for the decoupling of the substructure B from assembly AB:
A | |
\begin{bmatrix} Z | |
11 |
&
A | |
Z | |
12 |
&
A | |
0\\ Z | |
21 |
&
B | |
Z | |
22 |
&
B\\ 0& | |
Z | |
23 |
B | |
Z | |
32 |
&
B \end{bmatrix}- | |
Z | |
33 |
\begin{bmatrix} 0&0&0\\ 0&
B | |
Z | |
22 |
&
B | |
Z | |
23 |
\\ 0&
B | |
Z | |
32 |
&
B | |
Z | |
33 |
\end{bmatrix} =
A | |
\begin{bmatrix} Z | |
11 |
&
A | |
Z | |
12 |
&
A | |
0\\ Z | |
21 |
&
A | |
Z | |
22 |
&0\\ 0&0&0 \end{bmatrix}
Z\triangleq\begin{bmatrix} ZAB&0\\ 0&-ZB\end{bmatrix} ⇒ \tilde{Z}=LTZL
with:
L=\begin{bmatrix} I&0&0\\ 0&I&0\\ 0&0&I\\ 0&I&0\\ 0&0&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.
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} uAB=YABf+YABg\\ uB=-YBg \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} Y\triangleq&\begin{bmatrix} YAB&0\\ 0&-YB\end{bmatrix}\\[3pt] u=
AB | |
&\begin{bmatrix} u | |
1 |
AB | |
\\ u | |
2 |
AB | |
\\ u | |
3 |
B | |
\\ u | |
2 |
B | |
\\ u | |
3 |
\end{bmatrix};
AB | |
f\triangleq \begin{bmatrix} f | |
1 |
AB | |
\\ f | |
2 |
AB | |
\\ 0\\ 0\\ 0 \end{bmatrix}; g\triangleq \begin{bmatrix} 0\\ g | |
2 |
B | |
\\ 0\\ g | |
2 |
\\ 0 \end{bmatrix} \end{align}
In order to enforce compatibility, a similar approach is used as for the assembly task. Defining a
B
B=\begin{bmatrix} 0&-I&0&I&0 \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} u=Y(f-BT\boldsymbol{λ})\\ Bu=0 \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} coupling Y&\triangleq\begin{bmatrix}YA&0\ 0&YB\end{bmatrix} \\ decoupling Y&\triangleq\begin{bmatrix}YAB&0\ 0&-YB\end{bmatrix} \end{align}