r5.05.10 Dynamic analysis of structures with viscoelastic materials having frequency dependent properties#

Overview:

This reference document tackles the way how advanced capabilities are introduced in Code_Aster to take into account frequency dependent properties of viscoelastic materials. Computing of harmonic response and modes of vibration (real or complex) is addressed. Having frequency dependent modes is a step forward for the modal projection method and for model updating with an experimental modal basis as a reference. An iterative method is used in order to compute frequency dependent modes.

Finite element analysis#

State-of-the-art#

From a simulation point of view, the frequency response of a structure comprising at least one VEM is obtained by solving the dynamic equilibrium equation:

\(\left([{K}^{\text{*}}(\omega )]-{\omega}^{2}[M]\right)\left\lbrace u(\omega )\right\rbrace =\left\lbrace F(\omega )\right\rbrace\) , (2.1-1)

with \([{K}^{\text{*}}(\omega )]\) , the complex stiffness matrix as:

\([{K}^{\text{*}}(\omega )]=[{K}^{'}(\omega )]+i[{K}^{''}(\omega )]\) . (2.1-2)

The frequency dependence of the matrix comes from the use of frequency dependent complex moduli to represent viscoelastic behaviors [1] .

In many finite element codes, solving the system of Eq. () is not conventional because it needs to realize the stiffness matrix for each frequency step and to have computing procedures which are able to deal with frequency dependent matrices.

The direct response approach consists in solving Eq. () for each frequency step as:

\(\left\lbrace u(\omega )\right\rbrace ={\left([{K}^{\text{*}}(\omega )]-{\omega}^{2}[M]\right)}^{-1}\left\lbrace F(\omega )\right\rbrace \forall \omega \in [0,{\omega}_{max}]\) . (2.1-3)

This method has the advantage of computing the exact response of the system. But, as it is necessary to compute and inverse a complex matrix at each frequency step, the computing time can become prohibitive for industrial structures with several million degrees of freedom.

A few studies [2] [4] have shown the interest to solve Eq. () with dedicated modal response methods in order to improve the computational efficiency while maintaining the accuracy of the results. Modal response methods compute frequency responses by projecting the system of Eq. () on a modal basis \([T]\) , with the assumption:

\(\left\lbrace u(\omega )\right\rbrace =[T]\left\lbrace q(\omega )\right\rbrace\) . (2.1-4)

The projection of the model on the considered basis leads to a low order model, that decreases significantly the number of degrees of freedom and consequently the computing time of frequency responses:

\(\left\lbrace q(\omega )\right\rbrace ={\left({[T]}^{T}[{K}^{\text{*}}(\omega )][T]-{\omega}^{2}{[T]}^{T}[M][T]\right)}^{-1}{[T]}^{T}\left\lbrace F(\omega )\right\rbrace\) . (2.1-5)

Using properties of elastic models (real and frequency independent stiffness matrix), the standard basis of the popular spectral decomposition method combines normal modes solving of the classical eigenvalue problem:

\(\left([K]-{\omega}_{r}^{2}[M]\right)\left\lbrace {\Phi}_{r}\right\rbrace =\lbrace 0\rbrace\) (2.1-6)

and a static correction to ensure a correct representation of the low frequency contribution of truncated high frequency normal modes [5] . Classically, \([T]\) is composed by the normal modes between \(0\) and a minimum of \(1.5\times {\omega}_{max}\) .

However, frequency dependence of the VEM dynamic properties prevents from using the spectral decomposition method, because the modal basis \([\left\lbrace {\Phi}_{r=1,N}\right\rbrace ]\) is only valid at the frequency \({\omega}_{\mathit{ref}}\) , for which it has been computed. In practice, it has been demonstrated [3] the validity domain of the basis can be extended in a range around the frequency \({\omega}_{\mathit{ref}}\) . If the VEM of the model give increasing moduli with respect to the frequency, the use of \({\omega}_{\mathit{ref}}={\omega}_{max}\) will lead to a larger validity domain than for any other frequencies. Unfortunately, the validity domain of the basis may not extend over all the frequency range of interest. Inaccurate results may be obtained. To tackle this difficulty, static correction associated to the VEM can also be taken into account. This technique, presented in [U2.06.04], leads to very accurate results, but its implementation for an industrial case is tricky.

Computation of frequency dependent mode#

For a frequency dependent real stiffness matrix, \([K(\omega )]\) , the eigenvalue problem of Eq. () is written as:

\(\left([K(\omega )]-{\left({\omega}_{r}(\omega )\right)}^{2}[M]\right)\left\lbrace {\Phi}_{r}(\omega )\right\rbrace =\lbrace 0\rbrace\) , (2.2-1)

where the eigenfrequency, \({\omega}_{r}(\omega )\) , and the eigenvector, \(\left\lbrace {\Phi}_{r}(\omega )\right\rbrace\) , are frequency dependent too. For a fixed, given frequency, \(\omega ={\omega}_{p}\) , such eigensolutions can be obtained by solving the standard problem of Eq. (). Hence, the proposed method for solving the frequency dependent problem of Eq. () consists in using an iterative algorithm searching for:

\(∣{\omega}_{r}({\omega}_{p})-{\omega}_{p}∣<ϵ.{\omega}_{p}\forall {\omega}_{p}\in [{\omega}_{1,}{\omega}_{2}]\) , (2.2-2)

where \({\omega}_{1}\) and \({\omega}_{2}\) are respectively the lower and the upper cut-off frequencies of the eigenproblem and \(ϵ\) is the convergence criterion of the iterative algorithm. First, for \({\omega}_{p}={\omega}_{1}\) , the stiffness matrix is realized and the eigenproblem of Eq. () is solved. Then, \({\omega}_{p}\) is updated by taking the value of the \({n}^{\mathit{th}}\) eigenfrequency for which:

\(\mid {\omega}_{n}-{\omega}_{p}\mid =min\mid {\omega}_{r=1,N}-{\omega}_{p}\mid\) . (2.2-3)

Next, the stiffness matrix is realized for the new value of \({\omega}_{p}\) and the eigenproblem of Eq. () is solved again. The iterations will not stop while

\(∣{\omega}_{n}-{\omega}_{p}∣⩾ϵ.{\omega}_{p}\) . (2.2-4)

When the procedure converges, \({\omega}_{n}\) and \(\left\lbrace {\Phi}_{n}\right\rbrace\) are extracted to form the frequency dependent eigensolutions of Eq. () and the iterative algorithm will continue using \({\omega}_{p}={\omega}_{n+1}\) as the guess value to compute the next eigensolution. Finally, the algorithm will stop when \({\omega}_{p}>{\omega}_{2}\) . It means all the frequency dependent eigensolutions have been computed between \({\omega}_{1}\) and \({\omega}_{2}\) .

The computing time of the proposed method is classically driven by the number of frequency dependent modes to be computed, the number of degrees of freedom of the model and the value of the convergence criterion \(ϵ\) . It is also dependent on the number of normal eigenvalues which are computed as solutions of Eq. () at each iteration. Having all the eigenvalues between \({\omega}_{1}\) and \({\omega}_{2}\) for each iteration is not necessary and could lead to prohibitive computing times. In theory, a minimum of two eigenvalues may be sufficient to run iterations: \({\omega}_{n}\) to update \({\omega}_{p}\) as described in Eq. () and \({\omega}_{n+1}\) to continue when the procedure converges. In practice, this minimum number of eigenvalues will be determined by the capabilities of the numerical solver for Eq. (). It will be discussed for Code_Aster in section 3 .

The proposed method can be naturally extended to the study of damped structures comprising VEM. The frequency dependent eigenvalues and eigenvectors are then computed as solutions of the following problem:

\(\left([{K}^{\text{*}}({\lambda}_{r})]+{\lambda}_{r}^{2}[M]\right)\left\lbrace {\Psi}_{r}\right\rbrace =\left\lbrace 0\right\rbrace\) , (2.2-5)

with \({\lambda}_{r=1,N}\) the complex eigenvalues, and \(\left\lbrace {\Psi}_{r=1,N}\right\rbrace\) the associated complex modes. In this case, Eq. () is rewritten as:

\(∣\stackrel{̃}{\omega}{}_{r}({\omega}_{p})-{\omega}_{p}∣<ϵ.{\omega}_{p}\forall {\omega}_{p}\in [{\omega}_{1,}{\omega}_{2}]\) , (2.2-6)

with \(\stackrel{̃}{\omega}{}_{r}=∣{\lambda}_{r}∣\) , the corresponding eigenfrequency in a structural damping model [6] . \({\omega}_{p}\) will be updated by taking the value of the \({n}^{\mathit{th}}\) eigenfrequency for which:

\(\mid \stackrel{̃}{\omega}{}_{n}-{\omega}_{p}\mid =min\mid \stackrel{̃}{\omega}{}_{r=1,N}-{\omega}_{p}\mid\) . (2.2-7)

Computing frequency dependent real or complex modes should be a step forward for comparison with experimental modal bases. Indeed, when structures comprising VEM are tested, their frequency dependent behaviors are physically measured. The proposed method could help to validate or to update finite element models with experimental references.

Improvement of the modal projection method#

The frequency dependent modes can be used to form the basis of the modal projection method for response computations of VEM. They will extend the validity domain of the basis over all the frequency range of interest. In combination to a static correction, either real modes will be used for weakly damped structures, or complex modes for highly damped structures [3] . The static correction will be determined with the stiffness matrix realized for \(\omega =0\) . Real modes will be preferred to reduce computing times, since numerical solvers are much more efficient in this case. But for highly damped structures, the projection base composed with such modes may be insufficient to obtain accurate frequency responses. This can be improved using the modified Modal Strain Energy (MSE) method to reduce the errors [7] .

So, frequency dependent real modes, \(\left\lbrace \widehat{\Phi}{}_{r}(\omega )\right\rbrace\) , become solutions of a modified form of Eq. () in order to take into account the damping matrix (imaginary part of the stiffness matrix) as:

\(\left([{K}^{'}(\omega )]+\beta (\omega )[{K}^{''}(\omega )]-{\left(\widehat{\omega}{}_{r}(\omega )\right)}^{2}[M]\right)\left\lbrace \widehat{\Phi}{}_{r}(\omega )\right\rbrace =\lbrace 0\rbrace\) . (2.3-1)

\([{K}^{'}(\omega )]\) and \([{K}^{''}(\omega )]\) are defined by Eq. () and \(\beta (\omega )\) is calculated by the following empirical formula:

\(\beta (\omega )=\frac{\mathit{trace}[{K}^{''}(\omega )]}{\mathit{trace}[{K}^{'}(\omega )]}\) , (2.3-2)

where the trace for an \(N\times N\) matrix is defined as:

\(\mathit{trace}[A]=\sum_{j=1}^{N}{A}_{jj}\) . (2.3-3)

In the iterative algorithm, the use of the modified MSE method leads to calculate \(\beta ({\omega}_{p})\) for each iteration and to search for the solutions of the resulting eigenvalue problem:

\(\left([{K}^{'}({\omega}_{p})]+\beta ({\omega}_{p})[{K}^{''}({\omega}_{p})]-\widehat{\omega}{}_{r}^{2}[M]\right)\left\lbrace \widehat{\Phi}{}_{r}\right\rbrace =\lbrace 0\rbrace\) . (2.3-4)

Another method to improve the real modes could be the use of residual modes whose purpose is to represent the damping of VEM in the modal projection basis [8] . The coupling with frequency dependent modes has been presented in [U2.06.04].

Implementation in Code_Aster#

The approach to frequency dependent modes described above has been implemented in Code_Aster via a macro command, namely DYNA_VISCO. This script uses standard commands existing in Code_Aster for the conventional modal method. Python codes have been added for the new developments:

  • definition of a frequency dependent behavior,

  • computation of frequency dependent modes by the iterative algorithm,

  • frequency response computation with realization of the stiffness matrix at each frequency step.

In the iterations, modes are computed as solutions of Eq. () via the standard CALC_MODES command. The modal extraction method are described in [9] . Five eigensolutions are searched around \({\omega}_{p}\) for each iteration. This number is a good compromise between securing the convergence of the iterative algorithm and limiting the total computing time. In the case where the algorithm would not converge, the number of the searched eigensolutions will be automatically increased by the program.

The user can define several VEM with different behaviors and choose the frequency dependent modes to be used: real modes or beta-modes. At this time of the implementation, it is also possible to compute frequency dependent complex modes, but they cannot be used as a modal basis for response computation.

References#

    1. Nashif, D. I. G. Jones, J. P. Henderson, Vibration Damping , John Wiley and Sons (1985).

A-S. Plouin, E. Balmès, “Steel/Viscoelastic/Steel Sandwich Shells, Computational Methods and ExperimentalValidations”, IMAC (2000).

  1. Merlette, S. Germès, F. van Herpe, L. Jézéquel, D. Aubry, “The use of suitable modal bases for dynamic predictions of structures containing high damping materials”, ISMA , Leuven (2004).

  1. Abbadi, N. Merlette, N. Roy, “Response Computation of Structures with Viscoelastic Damping Materials using a Modal Approach - Description of the Method and Application on a Car Door Model Treated with High Damping Foam”, 5thSymposium on Automobile Comfort , Le Mans (2008).

  1. Géradin, D. Rixen, Théorie des Vibrations, Applications à la Dynamique des Structures , Editions Masson (1996).

D.J. Ewins, Modal Testing: Theory and Practice , Research Studies Press LTD. (1986).

B-G. Hu, M. A. Dokainish, W. M. Mansour, “A Modified MSE Method for Viscoelastic Systems: A Weighted Stiffness Matrix Approach”, ASME Journal of Vibration and Acoustics 117, 226-231 (1995).

  1. Bobillot, E. Balmès, “Iterative Techniques for Eigenvalue Solutions of Damped Structures Coupled with Fluids”, 43rd AIAA/ASME/ASCI/AHS/ASC Structures, Structural Dynamics and Materials Conference , Denver (2002).

  1. D’Haene, J. Lu, “Influence of PVB mechanical properties on the vibrational damping of a windscreen”, ISMA , Leuven (2002).