2. Finite element analysis#
2.1. 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\{u(\omega )\right\}=\left\{F(\omega )\right\}\), (2.1-1)
with \([{K}^{\text{*}}(\omega )]\), the complex stiffness matrix as:
\([{K}^{\text{*}}(\omega )]=[{K}^{\text{'}}(\omega )]+i[{K}^{\text{'}\text{'}}(\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\{u(\omega )\right\}={\left([{K}^{\text{*}}(\omega )]-{\omega }^{2}[M]\right)}^{-1}\left\{F(\omega )\right\}\forall \omega \in [\mathrm{0,}{\omega }_{\mathit{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 reverse 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\{u(\omega )\right\}=[T]\left\{q(\omega )\right\}\). (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 reduces the computing time of frequency responses:
\(\left\{q(\omega )\right\}={\left({[T]}^{T}[{K}^{\text{*}}(\omega )][T]-{\omega }^{2}{[T]}^{T}[M][T]\right)}^{-1}{[T]}^{T}\left\{F(\omega )\right\}\). (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\{{\Phi }_{r}\right\}=\{0\}\) (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 }_{\mathit{max}}\).
However, frequency dependence of the VEM dynamic properties prevents from using the spectral decomposition method, because the modal basis \([\left\{{\Phi }_{r=\mathrm{1,}N}\right\}]\) 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 gives increasing moduli with respect to the frequency, the use of \({\omega }_{\mathit{ref}}={\omega }_{\mathit{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. Consistently 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.
2.2. 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\{{\Phi }_{r}(\omega )\right\}=\{0\}\), (2.2-1)
Where the eigenfrequency, \({\omega }_{r}(\omega )\), and the eigenvector, \(\left\{{\Phi }_{r}(\omega )\right\}\), are frequency dependent too. For a fixed, given frequency, \(\omega ={\omega }_{p}\), such eigensolutions can be obtained by solving the standard problem of Eq. (). Consequently, 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}∣<ϵ\mathrm{.}{\omega }_{p}\forall {\omega }_{p}\in [{\omega }_{\mathrm{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 =\mathit{min}\mid {\omega }_{r=\mathrm{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}∣⩾ϵ\mathrm{.}{\omega }_{p}\). (2.2-4)
When the procedure converges, \({\omega }_{n}\) and \(\left\{{\Phi }_{n}\right\}\) 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\{{\Psi }_{r}\right\}=\left\{0\right\}\), (2.2-5)
with \({\lambda }_{r=\mathrm{1,}N}\) the complex eigenvalues, and \(\left\{{\Psi }_{r=\mathrm{1,}N}\right\}\) the associated complex modes. In this case, Eq. () is rewritten as:
\(∣\stackrel{̃}{\omega }{}_{r}({\omega }_{p})-{\omega }_{p}∣<ϵ\mathrm{.}{\omega }_{p}\forall {\omega }_{p}\in [{\omega }_{\mathrm{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 =\mathit{min}\mid \stackrel{̃}{\omega }{}_{r=\mathrm{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.
2.3. 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\{\widehat{\Phi }{}_{r}(\omega )\right\}\), 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}^{\text{'}}(\omega )]+\beta (\omega )[{K}^{\text{'}\text{'}}(\omega )]-{\left(\widehat{\omega }{}_{r}(\omega )\right)}^{2}[M]\right)\left\{\widehat{\Phi }{}_{r}(\omega )\right\}=\{0\}\). (2.3-1)
\([{K}^{\text{'}}(\omega )]\) and \([{K}^{\text{'}\text{'}}(\omega )]\) are defined by Eq. () and \(\beta (\omega )\) is calculated by the following empirical formula:
\(\beta (\omega )=\frac{\mathit{trace}[{K}^{\text{'}\text{'}}(\omega )]}{\mathit{trace}[{K}^{\text{'}}(\omega )]}\), (2.3-2)
Where the trace for an \(N\times N\) matrix is defined as:
\(\mathit{trace}[A]=\sum _{j=1}^{N}{A}_{\mathit{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}^{\text{'}}({\omega }_{p})]+\beta ({\omega }_{p})[{K}^{\text{'}\text{'}}({\omega }_{p})]-\widehat{\omega }{}_{r}^{2}[M]\right)\left\{\widehat{\Phi }{}_{r}\right\}=\{0\}\). (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].