3. Classical dynamic substructuring methods#
3.1. Introduction#
After studying the various stages of substructuring separately, and the techniques they involve, it seems interesting to present the main methods of dynamic substructuring: the Craig-Bampton method and that of Mac Neal and the method of reduced interface modes.
Craig-Bampton uses, as a basis for projecting substructures, constrained modes and normal modes with fixed interfaces [bib1].
In addition, MacNeal uses, as a basis for projecting substructures, attachment modes and normal modes with free interfaces [bib2].
We will discuss a method derived from these methods called reduced interface mode methods which uses substructures, interface modes and normal modes with blocked interfaces as a basis for projection.
3.2. Craig-Bampton method#
The following presentation only involves two substructures \({S}^{1}\) and \({S}^{2}\), but it can be generalized to any number of components. After studying each substructure separately, their projection bases (normal modes with fixed interfaces and constrained modes) are known. For each of them (identified by the exponent \(k\)), a partition of the degrees of freedom is established, distinguishing the vector of the internal degrees of freedom \({q}_{i}^{k}\) and the vector of the degrees of freedom of link \({q}_{j}^{k}\):
\({q}^{k}=\left\{\begin{array}{c}{q}_{i}^{k}\\ {q}_{j}^{k}\end{array}\right\}\)
Either:
\({\varphi }^{k}\) the matrix of the eigenvectors of the substructure \({S}^{k}\),
\({\psi }^{k}\) the matrix of the constrained modes of the \({S}^{k}\) substructure.
The projection base of \({S}^{k}\) is characterized by the matrix:
\({\Phi }^{k}=[{\varphi }^{k}{\psi }^{k}]\)
The transformation of RITZ (equation [eq]) allows us to write:
\({\eta }_{i}^{k}\) is the vector of the generalized degrees of freedom associated with the eigenmodes of \({S}^{k}\),
\({\eta }_{j}^{k}\) is the vector of the generalized degrees of freedom associated with the constrained modes of \({S}^{k}\).
However, normal modes are determined with fixed interfaces, and each constrained mode is obtained by imposing a unit displacement on a degree of link freedom, the others being blocked.
The generalized coordinates relating to the static deformations are then the values of the degrees of freedom of connection:
\({q}_{j}^{k}={\eta }_{j}^{k}\)
Let’s look at the contribution of substructure \({S}^{k}\) from an energy point of view. The kinetic and deformation energies are:
\({T}^{k}=\frac{1}{2}{q}^{{k}^{T}}\mathrm{.}{M}^{k}\mathrm{.}{\dot{q}}^{k}=\frac{1}{2}{({\Phi }^{k}{\dot{\eta }}^{k})}^{T}\mathrm{.}{M}^{k}\mathrm{.}{\Phi }^{k}{\dot{\eta }}^{k}=\frac{1}{2}{\dot{\eta }}^{k}\mathrm{.}{\stackrel{ˉ}{M}}^{k}\mathrm{.}{\dot{\eta }}^{k}\)
\({U}^{k}=\frac{1}{2}{q}^{{k}^{T}}\mathrm{.}{K}^{k}\mathrm{.}{\dot{q}}^{k}=\frac{1}{2}{({\Phi }^{k}{\dot{\eta }}^{k})}^{T}\mathrm{.}{K}^{k}\mathrm{.}{\Phi }^{k}{\dot{\eta }}^{k}=\frac{1}{2}{\dot{\eta }}^{k}\mathrm{.}{\stackrel{ˉ}{K}}^{k}\mathrm{.}{\dot{\eta }}^{k}\)
These expressions show the projections of the mass and stiffness matrices on the basis of the substructure. These matrices, called generalized matrices, verify a certain number of properties:
due to the orthogonality of the normal modes in relation to the stiffness and mass matrices, the upper left block of these matrices is diagonal. Subsequently, we will consider that these modes are normalized with respect to the mass matrix,
it can also be shown that the constrained modes are orthogonal to the normal modes with respect to the stiffness matrix [bib4].
The generalized stiffness and mass matrices therefore have the following form:
where \({\lambda }^{k}\) is the matrix of generalized rigidities associated with the eigenmodes of \({S}^{k}\).
The choice of the projection base with blocked interfaces therefore leads to a coupling of normal modes and static deformations by the mass matrix.
This method is interesting if a numerical study of substructures is envisaged. In fact, normal modes with fixed interfaces and constrained modes lend themselves well to calculation. On the other hand, their experimental determination is difficult, because it is difficult to achieve a good quality experimental embedment.
We also show in [bib 1] that this method is of order 2 in \(\omega /{\omega }_{m}\) where \({\omega }_{m}\) is the largest identified natural pulsation.
3.3. Mac Neal method#
It would be possible to present this method in a manner similar to that of Craig-Bampton, the only difference being in the use of normal modes with free interfaces and attachment modes. However, it seems interesting to adopt a slightly different approach, which makes it possible to arrive at a truncation criterion, and to make the attachment modes appear as the static contribution of the unidentified eigenmodes.
As before, the problem of a 2-component structure is considered. The method can be generalized to any number of substructures. In the following we will identify any quantity associated with the substructure \({S}^{k}\) by the exponent \(k\).
Here, we limit ourselves to the modal calculation of the global structure, so the external forces are zero, in order to alleviate theoretical developments. Also, the displacement of \({S}^{k}\) verifies the following dynamic equilibrium equation:
We will consider that the projection base of \({S}^{k}\) consists of all its modes specific to free interfaces. Thus, the dimension of the projected problem is equal to the dimension of the problem resulting from finite element modeling. In addition, we assume that a certain number of these modes have been determined, the others being unknown:
\({q}^{k}=[{\varphi }_{1}^{k}{\varphi }_{2}^{k}]\left\{\begin{array}{c}{\eta }_{1}^{k}\\ {\eta }_{2}^{k}\end{array}\right\}={\Phi }^{k}{\eta }^{k}\)
where:
\({\varphi }_{1}^{k}\) |
is the matrix of identified modal vectors from substructure \({S}^{k}\), |
\({\varphi }_{2}^{k}\) |
is the matrix of unidentified modal vectors from substructure \({S}^{k}\), |
\({\eta }_{1}^{k}\) |
is the vector of the generalized degrees of freedom associated with the identified eigenmodes of \({S}^{k}\), |
\({\eta }_{2}^{k}\) |
is the vector of the generalized degrees of freedom associated with the unidentified eigenmodes of \({S}^{k}\). |
The equation [eq] becomes, with the generalized coordinates defined above:
\(({\Phi }^{{k}^{T}}\mathrm{.}{K}^{k}\mathrm{.}{\Phi }^{k}-{\omega }^{2}{\Phi }^{{k}^{T}}\mathrm{.}{M}^{k}\mathrm{.}{\Phi }^{k}){\eta }^{k}={\Phi }^{{k}^{T}}{f}_{L}^{k}\)
either:
The natural modes are orthogonal to the mass and stiffness matrices and we choose to normalize them with respect to the mass matrix. So we have:
Let us now consider all of the two substructures. Each of them verifies the equations [eq] and [eq]. All of these dynamic equations constitute the following system:
By grouping the identified modes and the unidentified modes:
The equation for the transformation of RITZ becomes:
With these notations, the system of dynamic equations [eq] becomes:
This system of equations reflects the dynamic behavior of substructures separately. It does not represent the overall movements of the structure. To do this, it is necessary to add to it the conditions of connection between the two components.
The equations between substructures that link them derive from the equations [eq] and [eq], as well as from the organization of the base that we have chosen [eq]:
The equations [eq] and [eq] allow us to express the generalized coordinates relating to the unidentified modes:
From the equations [eq] and [eq], we therefore obtain:
Now, according to the formula [eq], we know that we can write bond matrices, in the form:
We therefore have, according to [eq] and [eq]:
We see the residual dynamic flexibility matrix associated with unidentified modes appear:
This term shows that the efforts to connect to the interface degrees of freedom are not zero if only known modes are taken into account in the modeling.
Therefore, the matrix problem [eq] can be reduced to the following equivalent system:
The unknowns of this problem are the generalized coordinates associated with the identified eigenmodes and the bond forces applied to the first substructure. The term residual flexibility, which also appears, is not known. A modeling that involves MacNeal’s attachment modes is proposed below.
There are two cases depending on whether or not the residual flexibility matrix is taken into account.
3.3.1. First case#
The residual dynamic flexibility matrix associated with unidentified modes is overlooked:
\({R}_{e}(\omega )=0\)
The resulting dynamic substructuring method is therefore very simple. It has the disadvantage of relying on a modal recombination method that is very sensitive to truncation effects.
3.3.2. Second case#
We use limited development: \({R}_{e}(\omega )={R}_{e}(0)+O({\omega }^{2}/{\omega }_{m}^{2})\).
Let’s adopt the following notations:
either:
\(n\) |
the number of modes in the complete database, |
\(m\) |
the number of unidentified modes, |
\(\Phi\) |
the matrix of \(n\) normal modes with free interfaces, |
\(\Psi\) |
the matrix of \(n\) attachment modes (defined for all degrees of freedom of the system, not only for the degrees of freedom of interface between the 2 substructures), |
\({\lambda }_{i}={\omega }_{i}^{2}\) |
the diagonal matrix of \(n\) eigenvalues. |
The complete matrix of attachment modes verifies: \(K\Psi =\mathrm{Id}\text{}\iff \text{}\Psi ={K}^{-1}\).
The complete modal base of normal modes with free interfaces constitutes an orthonormal basis. The stiffness matrix, expressed in this base, is written as:
\(\stackrel{ˉ}{K}={\Phi }^{T}\mathrm{.}K\mathrm{.}\Phi =\lambda\)
In the same way:
\({\stackrel{ˉ}{K}}^{-1}={\Phi }^{T}\mathrm{.}{K}^{-1}\mathrm{.}\Phi ={\lambda }^{-1}\)
A new expression for the complete matrix of attachment modes is therefore deduced:
\(\Psi ={K}^{-1}=\Phi \mathrm{.}{\lambda }^{-1}\mathrm{.}{\Phi }^{T}\)
Since the eigenmodes are orthogonal in pairs and the eigenvalue matrix is diagonal, the complete attachment mode matrix takes the final form:
\(\Psi =\sum _{i=1}^{n}{\varphi }_{i}{\lambda }_{i}^{-1}{\varphi }_{i}^{T}\)
Let us now consider the residual dynamic flexibility matrix, derived from the MacNeal method:
\({\mathrm{R}}_{e}(\omega )\mathrm{=}{\Phi }_{2}{({\lambda }_{2}\mathrm{-}{\omega }^{2}\mathrm{Id})}^{\mathrm{-}1}{\Phi }_{2}^{T}\mathrm{=}\mathrm{\sum }_{i\mathrm{=}m+1}^{n}{\varphi }_{i}{({\lambda }_{i}\mathrm{-}{\omega }^{2})}^{\mathrm{-}1}{\varphi }_{i}^{T}\)
When the number of modes identified is sufficiently large, the dynamic contribution becomes negligible against the static contribution:
\(\frac{{\omega }^{2}}{{\omega }_{i}^{2}}\ll 1\Rightarrow {R}_{e}(\omega )\approx {R}_{e}(0)=\sum _{i=m+1}^{n}{\varphi }_{i}{\lambda }_{i}^{-1}{\varphi }_{i}^{T}\)
where \({R}_{e}(0)\) is the residual static flexibility matrix.
To highlight the modal truncation effect, we can approach this matrix by expanding it to order 1:
\({R}_{e}(\omega )\approx \sum _{i=m+1}^{n}{\varphi }_{i}{\lambda }_{i}^{-1}(1+\frac{{\omega }^{2}}{{\omega }_{i}^{2}}){\varphi }_{i}^{T}\)
This matrix can be calculated according to the matrix of attachment modes:
\({R}_{e}(0)=\sum _{i=1}^{n}{\varphi }_{i}{\lambda }_{i}^{-1}{\varphi }_{i}^{T}-\sum _{i=1}^{m}{\varphi }_{i}{\lambda }_{i}^{-1}{\varphi }_{i}^{T}\Rightarrow {R}_{e}(0)=\Psi -\sum _{i=1}^{m}{\varphi }_{i}{\lambda }_{i}^{-1}{\varphi }_{i}^{T}\)
The second term in this formulation can be calculated exactly, since it only involves the modes (with free interfaces) identified.
Finally, note that in the MacNeal method, only the contribution of the attachment modes to the interface nodes is necessary.
The resolution of the [éq 3.3-4] system allows us to determine the natural frequencies of the global structure and the generalized coordinates of the natural modes. The way back to the expression of natural modes in the physical bases of substructures is achieved through the following relationship:
\({q}_{1}^{k}={\Phi }_{1}^{k}{\eta }_{1}^{k}+{R}_{e}^{k}(0){B}_{1}^{{k}^{T}}{f}_{{L}_{1}}^{T}\)
In the case of calculating the modes of the global structure, the MacNeal method therefore results in a problem with small eigenvalues. The mass and stiffness matrices resulting from the substructuring are symmetric. Two methods are in fact proposed, depending on whether residual flexibility is taken into account or not. The literature on the subject tends to show that the use of residual flexibility is essential to obtain reliable results [bib2].
3.4. Interface mode method (so-called reduced method)#
3.4.1. Introduction#
The interface mode method is based on the work of various authors; a complete theoretical justification can be found in [bib12]. We will not develop all these aspects in this document because the use made of them in Code_Aster is relatively simplified: it is just a matter of reducing the size of the reduced system by using the spectral properties of the dynamic operator condensed at the interface.
3.4.2. Define interface modes#
After studying each substructure separately to build the bases of normal modes, it is a question of considering the coupled problem to build the bases of interface modes. The procedure currently present in Code_Aster requires considering the problem coupled by dynamic substructuring. The Craig-Bampton method is therefore used without the normal modes with blocked interfaces, which is therefore equivalent to a static substructuration method by Guyan condensation. The following coupled system is obtained:
\({\stackrel{ˉ}{\mathrm{K}}}^{k}\mathrm{=}{\psi }^{{k}^{T}}\mathrm{.}{\mathrm{K}}^{k}\mathrm{.}{\psi }^{k}\text{}{\stackrel{ˉ}{\mathrm{M}}}^{k}\mathrm{=}{\psi }^{{k}^{T}}\mathrm{.}{\mathrm{M}}^{k}\mathrm{.}{\psi }^{k}\)
\(\psi\) refers to constrained modes here. The calculation of the eigenmodes of the global structure condensed at its interfaces by dynamic substructuring methods consists in solving a problem with the following matrix eigenvalues:
\((\stackrel{ˉ}{K}-{\omega }^{2}\stackrel{ˉ}{M}){\eta }_{R}+{L}^{T}\lambda =0\)
\(L{\eta }_{R}=0\)
It goes without saying that during this calculation of interface modes, it does not make sense to calculate as many modes as there are degrees of freedom of the interface because this would be equivalent to going back to the system of constrained modes. An error in truncating the projection base is therefore introduced at this level, but it remains physical and controllable.
The natural modes obtained are only defined on the degrees of interface freedom, it is therefore necessary to extend these modes to each of the coupled substructures. We thus obtain:
\({\psi }_{R}^{k}={\psi }^{k}{\eta }_{R}\)
3.4.3. Calculation of interface modes#
In practice, the calculation of interface modes with the method presented in paragraph 3.4.2 is not realistic. In fact, the operation of calculating constraint modes \(\psi\) is very expensive in memory and in terms of calculation time, especially since we then consider only a subfamily of \(N\) vectors, \(N\) being very low compared to the number of degrees of freedom of the interface. The approach adopted for the MODE_STATIQUE [U4.52.14] command is to build the interface modes using a preconditioner. If we note \(i\) the interface degrees of freedom, and \(c\) the complementary degrees of freedom, and the complementary degrees of freedom, by partitioning the stiffness matrix \(\mathrm{K}\), the constraint modes \(\psi\) satisfy
\(\left\{\begin{array}{c}{\psi }_{\mathrm{i}}\\ {\psi }_{\mathrm{c}}\end{array}\right\}\mathrm{=}\left\{\begin{array}{c}\mathit{Id}\\ \mathrm{-}{\mathrm{K}}_{\mathrm{cc}}^{\mathrm{-}1}{\mathrm{K}}_{\mathrm{ci}}\end{array}\right\}\)
Rather than building the stress modes directly on this principle, we assemble and build matrices of mass \(\tilde{{\mathrm{M}}_{\mathrm{ii}}}\) and stiffness \(\tilde{{\mathrm{K}}_{\mathrm{ii}}}\) defined only at the interface. These matrices are built on the basis of an Euler Bernoulli beam model using the connectivity of the interface nodes. Each pair of connected nodes of the interface is then connected by a cylindrical beam whose geometric properties are adjusted in relation to the shortest distance between two connected nodes of the interface. The diameter of the beam is then fixed at 20% of this minimum distance. The material used for the beams is steel. A mode family \(\tilde{{\phi }_{\mathrm{i}}}\) of this interface model is then calculated. Here, it is appropriate to calculate a number of \(M\) modes that is quite significantly greater than \(N\), where \(N\) is the number of interface modes that the user wishes to remember.
To calculate these modes, a first operation is carried out which makes it possible to determine the number of disjoint sub-parts of a given interface. For each disjoint part of the interface, the 6 rigid body modes are calculated. If we write \(d\) the number of disjoint parts, then the total number \(M\) in the calculated mode is equal to:
\(M\mathrm{=}3(E(\frac{N}{d})+2E(\frac{6}{d})+1)\) if \(E(\frac{N}{d})>7\)
\(M\mathrm{=}3(6+E(\frac{N}{d})+2)\)
where \(E(\mathrm{.})\) corresponds to the integer part of the quantity in question.
These \(M\) modes \(\tilde{{\phi }_{\mathrm{i}}}\) are then recorded statically on the model of the substructure, and then form the basis on which the interface modes are searched:
\(\tilde{\psi }\mathrm{=}\left\{\begin{array}{c}\tilde{{\psi }_{\mathrm{i}}}\\ \tilde{{\psi }_{\mathrm{c}}}\end{array}\right\}\mathrm{=}\left\{\begin{array}{c}\tilde{{\phi }_{i}}\\ \mathrm{-}{\mathrm{K}}_{\mathrm{cc}}^{\mathrm{-}1}{\mathrm{K}}_{\mathrm{ci}}\tilde{{\phi }_{i}}\end{array}\right\}\)
The mass and stiffness matrices of the macro element are then projected onto these vectors,
\(\stackrel{ˉ}{\mathrm{K}}\mathrm{=}{\tilde{\psi }}^{T}\mathrm{.}\mathrm{K}\mathrm{.}\tilde{\psi }\text{}\stackrel{ˉ}{\mathrm{M}}\mathrm{=}{\tilde{\psi }}^{T}\mathrm{.}\mathrm{M}\mathrm{.}\tilde{\psi }\)
and we are looking for their own modes:
\((\stackrel{ˉ}{\mathrm{K}}\mathrm{-}{\omega }^{2}\stackrel{ˉ}{\mathrm{M}}){\eta }_{R}\mathrm{=}0\)
The interface modes selected for the construction of the macro elements are finally the first \(N\) vectors defined by:
\(\psi \mathrm{=}\tilde{\psi }{\eta }_{\mathrm{R}}\)
The calculation of the matrices \(\tilde{{\mathrm{M}}_{\mathrm{ii}}}\) and \(\tilde{{\mathrm{K}}_{\mathrm{ii}}}\), and the associated modes \(\tilde{{\phi }_{\mathrm{i}}}\), is fast and inexpensive in memory, and makes it possible to limit the bearing operation to only the vectors of interest for constructing the interface modes. Of course, if we calculate as many interface modes \(\tilde{{\phi }_{\mathrm{i}}}\) as DDL, then the subspace described is identical to that described by the constraint modes.
3.4.4. Creation of reduced macro-elements#
For each of the substructures, we have the bases of normal modes and the bases of the interface modes. The creation of macro-elements is then identical to the case of the Craig-Bampton method. A partition of the degrees of freedom is established, distinguishing the vector of internal degrees of freedom \({q}_{i}^{k}\) and the vector of the degrees of freedom of link \({q}_{j}^{k}\):
\({q}^{k}=\left\{\begin{array}{c}{q}_{i}^{k}\\ {q}_{j}^{k}\end{array}\right\}\)
Either:
\({\phi }^{k}\) the matrix of the eigenvectors of the substructure \({S}^{k}\),
\({\psi }_{R}^{k}\) the interface mode matrix for the \({S}^{k}\) substructure.
The projection base of \({S}^{k}\) is characterized by the matrix:
\({\Phi }^{k}=[{\phi }^{k}{\psi }_{R}^{k}]\)
The transformation of RITZ (equation [eq]) allows us to write:
\({\eta }_{i}^{k}\) is the vector of the generalized degrees of freedom associated with the eigenmodes of \({S}^{k}\),
\({\eta }_{j}^{k}\) is the vector of the generalized degrees of freedom associated with the interface modes of \({S}^{k}\).
Let’s take a look at the contribution of component \({S}^{k}\) from an energy point of view. The kinetic and deformation energies are:
\({T}^{k}=\frac{1}{2}{q}^{{k}^{T}}\mathrm{.}{M}^{k}\mathrm{.}{\dot{q}}^{k}=\frac{1}{2}{({\Phi }^{k}{\dot{\eta }}^{k})}^{T}\mathrm{.}{M}^{k}\mathrm{.}{\Phi }^{k}{\dot{\eta }}^{k}=\frac{1}{2}{\dot{\eta }}^{k}\mathrm{.}{\stackrel{ˉ}{M}}^{k}\mathrm{.}{\dot{\eta }}^{k}\)
\({U}^{k}=\frac{1}{2}{q}^{{k}^{T}}\mathrm{.}{K}^{k}\mathrm{.}{\dot{q}}^{k}=\frac{1}{2}{({\Phi }^{k}{\dot{\eta }}^{k})}^{T}\mathrm{.}{K}^{k}\mathrm{.}{\Phi }^{k}{\dot{\eta }}^{k}=\frac{1}{2}{\dot{\eta }}^{k}\mathrm{.}{\stackrel{ˉ}{K}}^{k}\mathrm{.}{\dot{\eta }}^{k}\)
These expressions show the projections of the mass and stiffness matrices on the basis of the substructure. These matrices, called generalized matrices, verify a certain number of properties because of the orthogonality of the normal modes with respect to the stiffness and mass matrices; the upper left block of these matrices is diagonal. In addition, we will consider these modes to be normalized with respect to the mass matrix.
The generalized stiffness and mass matrices therefore have the following form:
where \({\lambda }^{k}\) is the matrix of generalized rigidities associated with the eigenmodes of \({S}^{k}\).
3.5. Implemented in Code_Aster#
3.5.1. Study of substructures separately#
The projection base of each substructure is composed of dynamic eigenmodes and static deformations.
The dynamic eigenmodes of the substructure are calculated with the classic Code_Aster operator: CALC_MODES [U4.52.02]. In the case of the Craig-Bampton substructure, the link interfaces should be blocked. This is done with the AFFE_CHAR_MECA [U4.44.01] operator.
The operator DEFI_INTERF_DYNA [U4.64.01] allows you to define the connection interfaces of the substructure. In particular, the type of the interface is specified, which can be either « CRAIGB » (Craig‑Bampton), or « MNEAL » (MacNeal), or finally « AUCUN » (no static modes calculated).
The operator DEFI_BASE_MODALE [U4.64.02] allows you to calculate the complete projection base of the substructure. Thus, the dynamic modes calculated previously are copied. Moreover, the static deformations are calculated according to the type defined in the operator DEFI_INTERF_DYNA [U4.64.01]. If the type is « CRAIGB « , the constrained modes of the interfaces of the substructure are calculated. If the type is « MNEAL « , the attachment modes of the substructure interfaces are calculated. If the type is « AUCUN « , no static deformation is calculated, which corresponds to a MacNeal base without static correction.
The operator MACR_ELEM_DYNA [U4.65.01] calculates the generalized stiffness and mass matrices of the substructure, as well as the bond matrices.
3.5.2. assembly#
The complete structure model is determined by the operator DEFI_MODELE_GENE [U4.65.02]. In particular, each substructure is defined by the macro-element that corresponds to it (from MACR_ELEM_DYNA) and the angles of rotation that allow it to be oriented. The links between substructures are defined by the data of the names of the two substructures involved and those of the two interfaces opposite each other. In the case of the interface mode method, it is necessary to specify the “REDUIT” option in order to allow coupling between macro-elements.
The numbering of the complete generalized problem is done by the operator NUME_DDL_GENE [U4.65.03]. The generalized mass and stiffness matrices of the complete structure are assembled according to this numbering with the operator ASSE_MATR_GENE [U4.65.04].