6. Transient response by classical dynamic substructuring#
6.1. Introduction#
The purpose of this reference documentation is to present the theoretical bases of the two methods for calculating transient response by dynamic substructuration available in Code_Aster. The first consists in carrying out a transitory calculation by substructuration for which the equations of the problem are projected onto the bases associated with each substructure. The difficulty lies in the double dualization of boundary conditions which results in a singular mass matrix. To use explicit integration schemes (which require the inversion of the mass matrix), it is therefore necessary to modify the treatment of the interfaces in the operator for calculating the transient response. The second method consists in determining the eigenmodes of the complete structure by substructuring and in projecting the equations of the transitory problem on this basis. The step of retrieving the final generalized results on a physical basis must therefore take into account this double projection.
The transient response calculation operator that receives the substructuring is the operator DYNA_TRAN_MODAL [U4.53.21]. Based on modal recombination methods, it was designed to solve transient problems in generalized coordinates and is very effective for large problems where it allows to reduce the number of degrees of freedom. On the other hand, it supports the taking into account of localized non-linearities (at the nodes) that one wishes to generalize in the case of substructuring.
In this report, we present the two methods of transient substructure calculation available in Code_Aster, as well as their implementation.
6.2. Transient calculation by projection on the basis of substructures#
6.2.1. Dynamic equations verified by the substructures separately#
Let’s say a \(S\) structure composed of \({N}_{S}\) substructures marked \({S}^{k}\). We assume that each substructure is modelled in finite elements. The vibratory behavior of substructures results from the external forces applied to them and from the bonding forces exerted on them by the other substructures. So for \({S}^{k}\), we have:
The unknown field of displacement, resulting from finite element modeling, is sought on an appropriate space, of reduced dimension (Ritz transformation) according to the expression:
The Ritz transformation [eq], applied to real and virtual displacement leads to the transient dynamic equation of the substructure [eq], makes it possible to write:
The problem defined by equation [eq] is symmetric. On the other hand, its dimension is determined by the number of modes taken into consideration (dynamic modes and static deformations). It is therefore necessary to solve a conventional transitory problem but of reduced size.
6.2.2. Dynamic equations verified by the global structure#
The matrix writing of the dynamic equation, verified by the global structure, is written simply from the dynamic equations verified by each substructure:
\(\begin{array}{}\left[\begin{array}{ccccc}{\stackrel{ˉ}{M}}^{1}& & & & \\ & \mathrm{...}& & & \\ & & {\stackrel{ˉ}{M}}^{k}& & \\ & & & \mathrm{...}& \\ & & & & {\stackrel{ˉ}{M}}^{{N}_{S}}\end{array}\right]\left\{\begin{array}{c}{\ddot{\eta }}^{1}\\ \mathrm{...}\\ {\ddot{\eta }}^{k}\\ \mathrm{...}\\ {\ddot{\eta }}^{{N}_{S}}\end{array}\right\}+\left[\begin{array}{ccccc}{\stackrel{ˉ}{C}}^{1}& & & & \\ & \mathrm{...}& & & \\ & & {\stackrel{ˉ}{C}}^{k}& & \\ & & & \mathrm{...}& \\ & & & & {\stackrel{ˉ}{C}}^{{N}_{S}}\end{array}\right]\left\{\begin{array}{c}{\dot{\eta }}^{1}\\ \mathrm{...}\\ {\dot{\eta }}^{k}\\ \mathrm{...}\\ {\dot{\eta }}^{{N}_{S}}\end{array}\right\}+\left[\begin{array}{ccccc}{\stackrel{ˉ}{K}}^{1}& & & & \\ & \mathrm{...}& & & \\ & & {\stackrel{ˉ}{K}}^{k}& & \\ & & & \mathrm{...}& \\ & & & & {\stackrel{ˉ}{K}}^{{N}_{S}}\end{array}\right]\left\{\begin{array}{c}{\eta }^{1}\\ \mathrm{...}\\ {\eta }^{k}\\ \mathrm{...}\\ {\eta }^{{N}_{S}}\end{array}\right\}\\ =\left\{\begin{array}{c}{\stackrel{ˉ}{f}}_{\mathrm{ext}}^{1}\\ \mathrm{...}\\ {\stackrel{ˉ}{f}}_{\mathrm{ext}}^{k}\\ \mathrm{...}\\ {\stackrel{ˉ}{f}}_{\mathrm{ext}}^{{N}_{S}}\end{array}\right\}+\left\{\begin{array}{c}{\stackrel{ˉ}{f}}_{L}^{1}\\ \mathrm{...}\\ {\stackrel{ˉ}{f}}_{L}^{k}\\ \mathrm{...}\\ {\stackrel{ˉ}{f}}_{L}^{{N}_{S}}\end{array}\right\}\end{array}\)
To which must be added the link equations:
\(\mathrm{\forall }k,l\text{}{\mathrm{L}}_{{S}^{k}\mathrm{\cap }{S}^{l}}^{k}\left\{{\eta }^{k}\right\}\mathrm{=}{\mathrm{L}}_{{S}^{k}\mathrm{\cap }{S}^{l}}^{l}\left\{{\eta }^{l}\right\}\text{}{\mathrm{f}}_{{L}_{{S}^{k}\mathrm{\cap }{S}^{L}}}^{k}\mathrm{=}\mathrm{-}{\mathrm{f}}_{{L}_{{S}^{k}\mathrm{\cap }{S}^{L}}}^{l}\)
This system can be written in the condensed form:
6.2.3. Double dualization of boundary conditions#
The condensed problem, given above, takes the form of a transient system to which a linear stress equation (in force and displacement) is associated. In Code_Aster, this type of problem is classical and it is solved by double dualization of boundary conditions [R3.03.01], i.e. by the introduction of auxiliary variables also called Lagrange multipliers to dualize boundary conditions. After the introduction of Lagrange multipliers, the matrix system takes the form:
where \({\lambda }_{1}\) and \({\lambda }_{2}\) are the Lagrange multipliers.
As we can see, the introduction of Lagrange multipliers makes the mass matrix singular. Therefore, the use of the explicit integration schemes developed in the operator DYNA_TRAN_MODAL [U4.53.21] of*Code_Aster* is impossible because they require the inversion of the mass matrix.
To make the matrix non-singular, it suffices to dualize with the same Lagrange multipliers, the condition on the second derivative of the link equations.
Thus the condition of continuity of movements [eq] is modified by the equivalent system:
where \(\eta °\) and \(\dot{\eta }°\) are the generalized movements and speeds at the initial moment.
The resulting matrix system is of the form [bib10]:
It can be seen that the mass matrix has the same shape as the stiffness matrix. They are therefore reversible. This system is therefore perfectly equivalent to the equation [eq] (it verifies the link conditions at any time) and it can be processed, in this form, by the operator DYNA_TRAN_MODAL.
6.2.4. Depreciation matrix treatment#
It can be seen that the condition of continuity of movements, formulated in equation [eq], results in an undamped second-order equation. When solving a transitory problem in a time step, any numerical error is likely to self-perpetuate, thus reducing the stability of the algorithm. To optimize the damping of the numerical error, it suffices to dualize the condition on the first derivative of the link equations with the same Lagrange multipliers multiplied by 2 (so as to make this damping critical):
The resulting matrix system is of the form [bib10]:
\(\left[\begin{array}{ccc}\mathrm{-}\text{Id}& \mathrm{L}& \text{Id}\\ {\mathrm{L}}^{\mathrm{T}}& \stackrel{ˉ}{\mathrm{M}}& {\mathrm{L}}^{\mathrm{T}}\\ \text{Id}& \mathrm{L}& \mathrm{-}\text{Id}\end{array}\right]\left\{\begin{array}{c}\ddot{{\lambda }_{1}}\\ \ddot{\eta }\\ \ddot{{\lambda }_{2}}\end{array}\right\}+\left[\begin{array}{ccc}\mathrm{-}2\text{Id}& 2\mathrm{L}& 2\text{Id}\\ 2{\mathrm{L}}^{\mathrm{T}}& \stackrel{ˉ}{\mathrm{C}}& 2{\mathrm{L}}^{\mathrm{T}}\\ 2\text{Id}& 2\mathrm{L}& \mathrm{-}2\text{Id}\end{array}\right]\left\{\begin{array}{c}\dot{{\lambda }_{1}}\\ \dot{\eta }\\ \dot{{\lambda }_{2}}\end{array}\right\}+\left[\begin{array}{ccc}\mathrm{-}\text{Id}& \mathrm{L}& \text{Id}\\ {\mathrm{L}}^{\mathrm{T}}& \stackrel{ˉ}{\mathrm{K}}& {\mathrm{L}}^{\mathrm{T}}\\ \text{Id}& \mathrm{L}& \mathrm{-}\text{Id}\end{array}\right]\left\{\begin{array}{c}{\lambda }_{1}\\ \eta \\ {\lambda }_{2}\end{array}\right\}\mathrm{=}\left\{\begin{array}{c}0\\ {\stackrel{ˉ}{\mathrm{f}}}_{\mathit{ext}}\\ 0\end{array}\right\}\)
It is therefore observed that the treatment of the numerical error on the link equations leads to the modification of the damping matrix of the transitory problem. This modification is completely comparable to that carried out on the mass matrix.
On the other hand, we did not want to generalize this treatment in the case of resolving unamortized transitional problems. This would have led us to create a temporary amortization matrix. One could have feared increasing calculation times, without any real benefit. Moreover, it is entirely possible for the user to define a damping matrix whose coefficients are zero. Since the modification of the latter is automatic, the non-damped transient system will be effectively solved, while optimizing the processing of any numerical error affecting the connection equations.
6.2.5. Treatment of initial conditions#
Consider a substructure \({S}^{k}\) characterized by its projection base \({\Phi }^{k}\) composed of normal modes and static deformations. It is assumed that initially the substructure \({S}^{k}\) is subject to a field of movement or speed (this does not modify the demonstration in any way) noted :math:``. The Ritz transformation allows us to write:
\({\mathrm{q}}_{0}^{k}\mathrm{=}{\Phi }^{k}{\eta }_{0}^{k}\)
where \({\eta }_{0}^{k}\) is the vector of the generalized displacements (or speeds) of \({S}^{k}\) to be determined.
The vector of the initial generalized displacements (or speeds) is determined as follows:
\({q}_{0}^{k}={\Phi }^{k}{\eta }_{0}^{k}\Rightarrow {\Phi }^{{k}^{T}}{q}_{0}^{k}=({\Phi }^{{k}^{T}}{\Phi }^{k})\cdot {\eta }_{0}^{k}\)
\(\mathrm{\Rightarrow }{\eta }_{0}^{k}\mathrm{=}{({\Phi }^{{k}^{T}}{\Phi }^{k})}^{\mathrm{-}1}{\Phi }^{{k}^{T}}{\mathrm{q}}_{0}^{k}\)
6.3. Transitional calculation on a global modal basis calculated by substructuring#
The second method developed consists in solving the transitory problem on the modal basis of the complete structure calculated by substructuring.
6.3.1. Calculation of the natural modes of the complete structure by substructuring#
Each substructure \({S}^{k}\) is represented by a projection base, composed of dynamic eigenmodes and static deformations, which we have noted \({\Phi }^{k}\). The projection base for the resulting complete structure is noted \(\Phi\).
The modal basis for the complete structure is calculated by substructuring. Each mode obtained is therefore a linear combination of the vectors of the substructure projection bases:
where:
\({\Phi }_{p}\) |
is the eigenmode matrix of the complete structure, |
\(\alpha\) |
is the generalized modal coordinate matrix of the structure. |
The projection of the matrices and vectors constituting the transitory problem, on the basis of the eigenmodes of the complete structure calculated by modal synthesis, makes it possible to determine:
|
\(\stackrel{ˉ}{\stackrel{ˉ}{M}}={\alpha }^{T}\stackrel{ˉ}{M}\alpha\) |
|
\(\stackrel{ˉ}{\stackrel{ˉ}{K}}={\alpha }^{T}\stackrel{ˉ}{C}\alpha\) |
|
\(\stackrel{ˉ}{\stackrel{ˉ}{C}}={\alpha }^{T}\stackrel{ˉ}{C}\alpha\) |
|
\({\stackrel{ˉ}{\stackrel{ˉ}{f}}}_{\mathrm{ext}}={\alpha }^{T}{\stackrel{ˉ}{f}}_{\mathrm{ext}}\) |
Because of the orthogonality of the natural modes of the structure calculated by modal synthesis, with respect to the matrices \(M\) and \(K\), the generalized mass and stiffness matrices obtained above are diagonal:
6.3.2. Dynamic equation verified by the global structure#
The complete structure is subject to external forces applied to it. So, we can write:
The unknown displacement field, resulting from finite element modeling, is replaced by its projection based on the natural modes of the structure, according to the formula:
where \(\eta\) is the vector of the generalized coordinates of the structure.
The Ritz transformation [eq], applied to the dynamic transient equation of the structure [eq], makes it possible to write:
The step of rendering on a physical basis requires taking into account the double projection (on the modal basis of the complete structure, then on the basis of projection of the substructures - cf. eq).
The problem defined by the equation [eq] is of a completely classical form. It is necessary to solve a symmetric transitory problem whose dimension is determined by the number of modes calculated by substructuring and whose mass and stiffness matrices are diagonal.
Finally, it should be noted that the treatment of the initial conditions is identical to the case of the transitory calculation by projection on the bases of the substructures (cf. § 6.2).
6.4. Comparative study of the two methods developed#
The theoretical bases, associated with the two methodologies implemented in Code_Aster to perform a transient response calculation using substructuring techniques, were presented in the previous chapters. Here, we specify their essential particularities.
The first methodology consists in calculating the transient response by substructuring. The equation verified by the complete structure is then projected onto the bases of the substructures. The precision of this method is therefore directly determined by the extent of these bases. The latter can be enriched without leading to prohibitive calculation times because the substructures are, in principle, relatively small in size. Be that as it may, it is difficult to estimate the modal truncation effect simply by knowing the modes of the bases of the substructures. On the other hand, the substructure projection bases are composed of modes that are not all orthogonal to each other (static natural and deformed modes). The generalized mass and stiffness matrices constituting the final transitory problem are therefore non-diagonal. Overall, their width of the band can be determined from the number of static deformations of the projection bases of the substructures. The duration of integration in the transient calculation operator DYNA_TRAN_MODAL will therefore be all the longer as there are degrees of interface freedom. On the other hand, the maximum integration time step admissible by the explicit integration schemes is determined from the maximum frequency of the projection base. In the case of a transient calculation by substructuration, this frequency results, in principle, from static modes whose diagonal terms are high in the generalized stiffness matrix and low in the generalized mass matrix. Therefore, the integration time step cannot be determined a priory. Experience shows that it is very low, with respect to the natural frequencies of the bases of the substructures and that the use of the DYNA_TRAN_MODAL adaptive time step integration scheme is very profitable.
The second methodology consists in making a transitory calculation on the modal basis of the complete structure obtained by substructuring. It is known that the step consisting in calculating the natural modes of the structure can be expensive in terms of calculation time. This is all the more true when considering non-linear forces because the projection base must then be sufficiently extended to properly represent the dynamics of the system. On the other hand, the modal base on which the transitory response is calculated is generally of a smaller dimension than that determined by the eigenvectors of the substructures (eigenmodes and static deformations). Double projection therefore amounts to introducing a cutoff frequency. It is therefore to be expected that this method will be less accurate than the previous one. However, the calculation of natural modes makes it possible to estimate the modal truncation effect. On the other hand, it can be used to validate substructure models if experimental results are available. Finally, the essential advantage of this method is that the mass and stiffness matrices used in the transitory calculation are diagonal. Digital integration is therefore very fast.
6.5. Implemented in Code_Aster#
6.5.1. Study of substructures separately#
If you want to introduce Rayleigh damping, the parameters
and
of this damping are defined by the operator DEFI_MATERIAU [U4.43.01].
The treatments of the substructures are identical to the case of modal and harmonic calculation. Dynamic eigenmodes are calculated with the operator CALC_MODES [U4.52.02]. The conditions at the link interfaces are applied with the operator AFFE_CHAR_MECA [U4.44.01].
The operator DEFI_INTERF_DYNA [U4.64.01] allows you to define the connection interfaces of the substructure. The operator DEFI_BASE_MODALE [U4.64.02] makes it possible to calculate the complete projection base of the substructure (copying the natural modes and calculating the static deformations).
The operator MACR_ELEM_DYNA [U4.65.01] calculates the generalized stiffness, mass, and possibly damping matrices of the substructure, as well as the connection matrices. Rayleigh damping is taken into account by completing operand MATR_AMOR. Proportional damping is introduced by operand AMOR_REDUIT.
The transitory load is defined, at the level of the substructure, by the operators AFFE_CHAR_MECA [U4.44.01] (application of the force on the mesh), CALC_VECT_ELEM [U4.61.02] (calculation of the associated elementary vectors) and ASSE_VECTEUR [U4.61.23] (assembly of the load vector on the mesh of the substructure).
The CREA_CHAMP [U4.72.04] operator, which allows you to assign a field to the nodes of a model, allows you to describe the initial displacement field and/or the initial speed field of the substructure.
6.5.2. Generalized model assembly#
As in the case of modal and harmonic calculation, the complete structure model is defined by the operator DEFI_MODELE_GENE [U4.65.02]. Its numbering is done by the operator NUME_DDL_GENE [U4.65.03]. The generalized mass, stiffness and possibly damping matrices of the complete structure are assembled according to this numbering with the operator ASSE_MATR_GENE [U4.65.04].
The loads are projected onto the bases of the substructures to which they are applied, then assembled using the numbering from NUME_DDL_GENE [U4.65.03] by the operator ASSE_VECT_GENE [U4.65.05].
The initial generalized displacements and the initial generalized speeds for each substructure are calculated by the operator ASSE_VECT_GENE [U4.65.05]. This operator also assembles these vectors according to the numbering from NUME_DDL_GENE [U4.65.03].
In the case of a transitory calculation projected onto the « bases » of the substructures, the generalized matrices and assembled vectors obtained at the end of this step are directly used for the transitory calculation. In the case of a calculation on the modal basis of the complete structure calculated by substructuring, specific operations must be carried out which are presented in § 5.3.
6.5.3. Calculation of the modal base of the complete structure and projection#
This chapter is specific to the transitory calculation on a modal basis calculated by substructuring.
The modal base of the complete structure is calculated with the classical Code_Aster operator: CALC_MODES [U4.52.02]. A numbering of the final generalized problem is defined with the operator NUME_DDL_GENE [U4.65.03]. The generalized mass, stiffness and possibly damping matrices are projected based on the natural modes of the structure with the operator PROJ_MATR_BASE [U4.63.12]. The generalized vectors corresponding to external loads are projected based on the eigenmodes of the structure with the operator PROJ_VECT_BASE [U4.63.13].
6.5.4. Resolution and retrieval on a physical basis#
The calculation of the transient response of the complete structure is performed by the operator DYNA_TRAN_MODAL [U4.53.21].
The return of results on a physical basis involves the operator REST_GENE_PHYS [U4.63.31]; it is identical to the case of modal and harmonic calculation. You can use the DEFI_SQUELETTE [U4.24.01] operator to create a « skeleton » mesh. Coarser than computational meshing, it makes it possible to reduce the time required for graphics processing.
6.6. Conclusion#
In this report, we presented the work carried out to introduce, in*Code_Aster*, the calculation of linear transient response by dynamic substructuring. The methods that have been chosen consist, for the first of them, in projecting the transient equations onto the « bases » of each substructure, composed of dynamic eigenmodes and static deformations, and for the second, in calculating the eigenmodes of the complete structure by substructuring and in projecting the transient equations there.
We started with an overview of the theoretical bases on which the first method of transitory substructuration is based to arrive at the matrix formulation of the final problem. In particular, an original treatment of the displacement continuity equation to make the mass matrix invertible and to ensure optimal stability of the integration algorithm, led us to modify the shape of the mass and damping matrices of the transient problem.
For the second method, the essential difficulty consists in restoring the results obtained in generalized coordinates on a physical basis. Indeed, it is necessary to take into account the double projection: on the modal basis of the complete structure on the one hand, and on the basis of substructures on the other hand.
The developments carried out resulted in changes to the operators DYNA_TRAN_MODAL [U4.53.21] and REST_GENE_PHYS [U4.63.31]. Their syntax has been modified very little, so that their use is identical during a substructure calculation and a direct calculation by modal recombination.