# Vibrations non linéaires forcées de plaques minces amorties par une méthode asymptotique numérique

Faiza Boumediene, Jean-Marc Cadou, Abdelhamid Miloudi, Laëtitia Duigou, El Hassan Boutyour

► To cite this version:

Faiza Boumediene, Jean-Marc Cadou, Abdelhamid Miloudi, Laëtitia Duigou, El Hassan Boutyour. Vibrations non linéaires forcées de plaques minces amorties par une méthode asymptotique numérique. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01416068

HAL Id: hal-01416068

https://hal.science/hal-01416068v1

Submitted on 14 Dec 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Vibrations non linéaires forcées de plaques minces amorties par une méthode asymptotique numérique

F. Boumediene<sup>a</sup>, J. M. Cadou<sup>b</sup>, A. Miloudi<sup>a</sup>, L. Duigou<sup>b</sup>, E. H. Boutyour<sup>c</sup>

<sup>a</sup> Laboratoire de Mécanique Avancée, Faculté de Génie Mécanique & Génie des Procédés, USTHB, BP 32, El Alia, 16111 Bab Ezzouar, Alger, Algérie.
boumedienefaiza@yahoo.fr; amiloudi@yahoo.fr.

<sup>b</sup> Laboratoire d'Ingénierie des Matériaux de Bretagne, Université Européenne de Bretagne, Université de Bretagne Sud, Rue de Saint Maudé, BP 92116, 56321 Lorient cedex, France.
jean-marc.cadou@univ-ubs.fr; laetitia.duigou@univ-ubs.fr.

<sup>c</sup> Département de Physique Appliquée, Faculté des Sciences et Techniques, Université Hassan I, BP 577, Settat, Maroc.
boutyour@lpmm.sciences.univ-metz.fr.

**Résumé** Ce travail concerne la résolution de problèmes de vibrations forcées non linéaires de plaques rectangulaires amorties élastiques minces soumises à une excitation harmonique par une méthode asymptotique numérique. Cette dernière est une méthode de continuation qui permet de transformer le problème non linéaire en une séquence de problèmes linéaires ayant tous la même matrice de rigidité. Donc, une seule triangulation de matrice permet d'obtenir une grande partie de la solution. Des résultats numériques sont rapportés et comparés à ceux de la littérature.

**Mots-clés** Vibration non linéaire, amortissement, méthode asymptotique numérique, équilibrage harmonique, approximation de Padé, plaques.

## 1 Introduction

Pour réduire le poids des structures, les plaques minces sont largement répandus dans beaucoup de domaines industriels tels que l'aéronautique, la mécanique, le génie civil, etc... Ces structures sont souvent forcées à vibrer avec de grandes amplitudes induisant des non linéarités significatives de type géométriques. Ces dernières sont des sources de phénomènes d'instabilité complexes. Pour réduire ces effets, une approche commune consiste à introduire un amortissement dans la structure. La difficulté principale dans l'analyse de tels problèmes se situe à la présence, d'une part, de la non linéarité qui peut conduire à des phénomènes complexes telles que les bifurcations et dans d'autre part à l'amortissement qui mène généralement à des solutions complexes.

Dans la littérature, les principales études réalisées tenant compte de l'amortissement sont généralement limitées à la vibration linéaire conduisant alors à la résolution d'un problème de valeurs propres complexes. Seulement quelques travaux associent une non linéarité géométrique et l'amortissement: Amabili a étudié expérimentalement et numériquement des vibration des plaques avec un amortissement visqueux et soumises à une excitation harmonique [2] et [4] et numériquement des panneaux cylindrique circulaires avec un amortissement visqueux [3]. Ganapathi et al ont déterminé le facteur de perte d'une poutre sandwich avec une couche centrale viscoélastique; en utilisant la méthode des éléments finis et une approche itérative [13]. Daya et al [11] ont proposé une équation en fonction de l'amplitude de déplacement basée sur
une méthode d'équilibrage harmonique associée à la méthode de Galerkin pour étudier les vibrations des poutres et des plaques sandwichs avec des couches centrales viscoélastiques. La même approche a été prolongée à l'analyse d'anneaux circulaires sandwichs par Boutyour et al [8]. Touzé et al ont établi des modèles réduits pour l'étude de systèmes géométriquement non linéaires amortis [25]. Ils ont considéré un amortissement visqueux modal. Ribeiro et al [19] et [21] ont employé la méthode des éléments finis hiérarchique pour étudier profondément la réponse non linéaire de plaques rectangulaires complètement encastrées. Dans ces derniers articles, les auteurs ont considéré l'amortissement de type hystérésis proportionnel à la masse. Une approche semblable a été employée pour étudier la réponse forcée des plaques simplement appuyées avec les bords immobiles par Ribeiro [24] ainsi que des poutres et des plaques (minces et épaisses) Ribeiro et al [23]. Dans ces derniers travaux, les auteurs ont considéré un amortissement proportionnel à la rigidité de la structure.

Le but de ce travail est de traiter les vibrations non linéaires de plaques amorties par une méthode asymptotique numérique. En utilisant la théorie de plaque de Von Karman, la méthode d'équilibrage harmonique et le principe de Hamilton, l'équation de mouvement initiale est transformée en une forme adaptée à la méthode asymptotique. En couplant la technique de perturbation et la méthode des éléments finis, le problème non linéaire est transformé en une séquence de problèmes linéaires algébriques. Afin d'augmenter la longueur de pas et de réduire ainsi le temps de calcul, la méthode des approximants de Padé est utilisée [12]. Les résultats obtenus sont comparés à la littérature et aux approches analytiques. Pour montrer l'efficacité de l'algorithme proposé, plusieurs exemples de plaques rectangulaires avec diverses excitations et conditions aux limites sont présentés.

## 2 Formulation Mathématique

Considérons une plaque rectangulaire mince avec un système de coordonnée (o, x, y, z), l'origine est située en un coin de la plaque. Les composantes de déplacement de la surface moyenne de la plaque sont notées u, v et w, où u et v sont les déplacements de membrane et w le déplacement transversal dans les directions x, y et z respectivement. La déformation associée de Green-Lagrange varie linéairement avec l'épaisseur z:

$$ \epsilon = \gamma + zk \quad (1) $$

où $\gamma = \gamma^1 + \gamma^{nl}$ est la déformation généralisée de membrane qui peut être décomposée en une partie linéaire et une partie non linéaire et k est la déformation de flexion.

Les déformations sont liées aux déplacements suivant la théorie de plaque de Von-Karman comme suit:

$$ \gamma^1 = \left\{ \begin{matrix} \partial u / \partial x \\ \partial v / \partial y \\ \partial u / \partial y + \partial v / \partial x \end{matrix} \right\} , \quad \gamma^{nl} = \frac{1}{2} \left\{ \begin{matrix} (\partial w / \partial x)^2 \\ (\partial w / \partial y)^2 \\ 2(\partial w / \partial x)(\partial w / \partial x) \end{matrix} \right\} , \quad k = - \left\{ \begin{matrix} \partial^2 w / \partial x^2 \\ \partial^2 w / \partial y^2 \\ 2\partial^2 w / \partial x \partial y \end{matrix} \right\} \quad (2) $$

Pour des raisons de simplicité de lecture, la déformation est écrite sous la forme suivante:

$$ \epsilon = \left[ B_1 + \frac{1}{2} B_{nl}(U) \right] U \quad (3) $$

où U est le vecteur de déplacement donné par $U = \{u, v, w\}^t$

Les contraintes généralisées N sont liées à la déformation, pour un matériau homogène et isotrope, par la relation:

$$ N = \left\{ \begin{matrix} N_x \\ N_y \\ N_{xy} \end{matrix} \right\} = D \gamma \quad (4) $$
où D est une matrice symétrique contenant les coefficients élastiques du matériau.

En négligeant les moments d'inertie, l'énergie cinétique est donnée par:

$$T = \frac{1}{2} \int_{S_0} \rho h \dot{U}^2 dS \qquad (5)$$

où $S_0$ est la surface moyenne de la plaque, $\rho$ est la masse volumique, $h$ est l'épaisseur de la plaque et le point signifie une dérivation par rapport au temps $t$.

En négligeant la contrainte transversal $\sigma_z$ sous les hypothèses de Kirchhoff, l'énergie potentielle $V$ d'une plaque est donnée par:

$$V = \frac{1}{2} \int_{S_0} N \varepsilon \, dS \qquad (6)$$

Les expressions de (1) à (4) sont alors introduites dans l'équation (6). En Tenant compte des coefficients d'amortissement et en utilisant le principe de Hamilton, on obtient l'équation de mouvement sous la forme suivante:

$$\left\{ \begin{array}{l} \displaystyle \int_{S_0}^{t} B(U) \cdot N \, dS + C \dot{U} - M \ddot{U} = F \\[1em] \displaystyle N = D \left[ B_l + \frac{1}{2} B_{nl} (U) \right] U \end{array} \right. \qquad (7)$$

où $M$ est la matrice de masse, $F$ est le vecteur d'excitation extérieure et $C$ est la matrice d'amortissement visqueux de type Rayleigh $C = \alpha M + \beta K$ ($\alpha$ et $\beta$ sont deux paramètres et $K$ est la matrice de rigidité élastique linéaire).

$$\text{et } B(U) = B_l + B_{nl}(U)$$

Le problème est de résoudre le système (7) où les inconnus sont les déplacements, les contraintes et les pulsations.

### 3 Méthode de l'équilibrage harmonique

L'excitation harmonique considérée est donnée par:

$$F(t) = \sum_{j=0}^{H-1} F_c^j \cos j\omega t + F_s^j \sin j\omega t \qquad (8)$$

où $H$ est le nombre d'harmoniques.

On suppose alors que la réponse de la plaque est harmonique et peut être écrite sous la forme suivante:

$$\left\{ \begin{array}{l} U(t) = \sum_{j=0}^{H-1} U_c^j \cos j\omega t + U_s^j \sin j\omega t \\ N(t) = \sum_{j=0}^{H-1} N_c^j \cos j\omega t + N_s^j \sin j\omega t \end{array} \right. \qquad (9)$$

On introduit de nouveaux vecteurs $\bar{U}$ et $\bar{N}$ contenant toutes les harmoniques de la relation (9) et définis par:

$$\left\{ \begin{array}{l} \bar{U} = [U_c^0, U_c^1, U_s^1, \dots, U_c^i, U_s^i, \dots, U_c^{H-1}, U_s^{H-1}] \\ \bar{N} = [N_c^0, N_c^1, N_s^1, \dots, N_c^i, N_s^i, \dots, N_c^{H-1}, N_s^{H-1}] \end{array} \right. \qquad (10)$$

où "c" note le facteur de co-sinus et "s" le sinus.

En insérant la série (9) dans le système (7) et en utilisant la méthode de l'équilibrage harmonique, le problème est alors réécrit comme suit:
$$
\left\{
\begin{array}{l}
\displaystyle \int_{S_0}^t \overline{B}(\overline{U}) \overline{N} \, dS + \omega \cdot \overline{C} \overline{U} - \omega^2 \overline{M} \overline{U} = \overline{F} \\[1em]
\displaystyle \overline{N} = \overline{D} \left[ \overline{B}_1 + \frac{1}{2} \overline{B}_{nl} (\overline{U}) \right] \overline{U}
\end{array}
\right.
\qquad (11)
$$

où les matrices $\overline{M}$, $\overline{C}$, $\overline{D}$, $\overline{B}^1$ et $\overline{B}^{nl}$ sont calculées à partir des matrices M, C, D, F, $B^1$ et $B^{nl}$ respectivement (cf. l'annexe). Le vecteur $\overline{F}$ est calculé à partir du vecteur F et est écrit sous la forme (10).

Ce dernier système d'équation (11) est non linéaire. Notre objectif est de le résoudre par une
méthode asymptotique numérique. Pour cela, nous l'avons écrit sous une forme quadratique par
rapport à un vecteur inconnu mixte $\Lambda=^t[\overline{U}, \overline{N}]$ et des pulsations:

$$
\langle L\Lambda, \delta\Lambda \rangle - \omega^2 \langle \overline{M}\Lambda, \delta\Lambda \rangle + \Omega \langle \overline{C}\Lambda, \delta\Lambda \rangle + \langle Q(\Lambda, \Lambda), \delta\Lambda \rangle = \langle \overline{F}, \delta\Lambda \rangle \quad (12)
$$

où L(.) est un opérateur linéaire et Q(..) est un opérateur quadratique définis par les
expressions:

$$
\langle L\Lambda, \delta\Lambda \rangle = \int_{S_0}^{t} \overline{B}_1 \cdot \overline{N} \, dS
$$

$$
\langle Q(\Lambda, \Lambda), \delta\Lambda \rangle = \int_{S_0}^{t} \overline{B}_{nl}(\overline{U}) \overline{N} \, dS
$$

$\Omega = \omega^2$

4 Méthode Asymptotique Numérique (MAN)

4.1 Développement asymptotique

Considérons une solution régulière $(\Lambda_0, \omega_0, \Omega_0)$ du problème non linéaire (12). L'idée fondamentale de la MAN consiste à chercher la solution à proximité de ce point sous la forme d'une série entière par rapport à un paramètre de chemin 'a':

$$
\left\{
\begin{aligned}
\Lambda(a) &= \Lambda_0 + a\Lambda_1 + a^2\Lambda_2 + \cdots + a^n\Lambda_n \\
\omega(a) &= \omega_0 + a\omega_1 + a^2\omega_2 + \cdots + a^n\omega_n \\
\Omega(a) &= \Omega_0 + a\Omega_1 + a^2\Omega_2 + \cdots + a^n\Omega_n
\end{aligned}
\right.
\qquad (13)
$$

où $(\Lambda_p, \omega_p, \Omega_p)$ représente l'inconnue à l'ordre de troncature p du développement asymptotique (13).

Le paramètre de chemin utilisé dans la série (13) peut être identifié comme étant la projection de l'incrément de déplacement $(\overline{U} - \overline{U}_0)$ et de l'incrément de fréquence $(\omega - \omega_0)$ sur le vecteur tangent $(\overline{U}_1, \omega_1)$:

$$
a = \langle \overline{U} - \overline{U}_0, \overline{U}_1 \rangle + (\omega - \omega_0) \omega_1 \quad (14)
$$

où <.,.> désigne le produit scalaire euclidien.

En introduisant les équations (13) dans l'équation (12) et (14) et en identifiant suivant les
puissances de 'a', on obtient un ensemble de problèmes linéaires:

ordre 1

$$
\left\{
\begin{array}{l}
L_t(\Lambda_1) = (\Omega_1 \overline{M} - \omega_1 \overline{C}) \Lambda_0 \\
\int^t \overline{U}_1 \overline{U}_1 + \omega_1^2 = 1
\end{array}
\right.
\qquad (15)
$$
ordre p

$$
\left\{
\begin{array}{l}
L_t \Lambda_p = \sum_{i=0}^{p-1} (\Omega_{p-i} \overline{M} - \omega_{p-i} \overline{C}) \Lambda_i - \sum_{i=1}^{p-1} Q(\Lambda_i, \Lambda_{p-i}) \\
\\
\quad + \omega_1 \omega_p = 0
\end{array}
\right.
\qquad (16)
$$

où :

$$
L_t(\cdot) = L(\cdot) + 2.Q(\Lambda_0, \cdot) + \omega_0 \bar{C} - \omega_0^2 \bar{M}
$$

4.2 Discrétisation par la méthode des éléments finis

Après linéarisation du problème, on revient à l'écriture en déplacement pure pour la
résolution avec la méthode d'éléments finis classique en déplacement. Après discrétisation en
éléments finis, on obtient les étapes de résolution comme suit:

Ordre 1

$$
\left\{
\begin{array}{l}
[\bar{K}_t]\{q_1\} = \omega_1[2\omega_0\bar{M} - \bar{C}]\{q_0\} \\
\langle q_1, q_1 \rangle + \omega_1^2 = 1 \\
\bar{N}_1 = \bar{D} \bar{B}(q_0)q_1 \\
\Omega_t = 2\omega_0\omega_1
\end{array}
\right.
\qquad (17)
$$

Ordre p

$$
\left\{
\begin{aligned}
[\bar{K}_t]\{q_p\} &= \omega_p[2\omega_0\bar{M} - \bar{C}]\{q_0\} + \bar{F}_p^{nl} \\
\langle q_1, q_p \rangle + \omega_1\omega_p &= 0 \\
\bar{N}_p &= \bar{D}\bar{B}(q_0)q_p + \frac{1}{2}\bar{D}\sum_{i=1}^{p-i} \bar{B}^{nl}(q_{p-i})q_i \\
\Omega_p &= 2\omega_0\omega_p + \underbrace{\sum_{i=1}^{p-1} \omega_i\omega_{p-i}}_{\Omega_p^{nl}}
\end{aligned}
\right.
\qquad (18)
$$

où $[\bar{K}_t]$ est la matrice tangente au point de départ $(\Lambda_0, \omega_0, \Omega_0)$ et $\{q_p\}$ est la forme discrétisée du déplacement $\bar{U}_p$

$$
F_p^{nl} = \Omega_p^{nl} \cdot \overline{M} \cdot q_0 + \sum_{i=1}^{p-1} (\Omega_{p-i} \overline{M} - \omega_{p-i} \overline{C}) \cdot q_i - \int_{S_0} \left[ \overline{B}(q_0) \overline{N}_p^{nl} + \sum_{i=1}^{p-1} (\overline{B}_{nl}(q_i) \overline{N}_{p-i}) \right] dS
$$

Les plaques sont modélisées avec des éléments coques triangulaires classique DKT [6] qui ont
trois nœuds et cinq degrés de libertés par nœud (u, v, w, θx, θy).

Tous les paramètres inconnus des séries (13) peuvent être déterminés en résolvant successivement les équations (17) et (18), en chaque ordre p.

4.3 L'intervalle de validité et les approximants de padé

Les solutions polynomiales (13) coïncident presque parfaitement en dessous du rayon de
convergence mais elles divergent hors de cette zone de validité. Cette limite peut être calculée
automatiquement en utilisant un critère simple proposé par Cochelin [10]:
$$ a = \left( \eta \frac{\|q_1\|}{\|q_n\|} \right)^{\frac{1}{n-1}} \qquad (19) $$

où $η$ est un petit nombre donné.

La représentation polynomiale (13) peut être améliorée en utilisant les fractions rationnelles appelées les approximants de Padé [12]:

$$ \left\{ \begin{aligned} q(a) &= q_0 + \sum_{i=1}^{n-1} f_i(a) a^i q_i \\ \omega(a) &= \omega_0 + \sum_{i=1}^{n-1} f_i(a) a^i \omega_i \end{aligned} \right. \qquad (20) $$

où $f_i(a)$ sont les fractions rationnelles ayant le même dénominateur.

La zone de validité de l'approximation rationnelle (20) est déterminée en considérant la différence relative entre les déplacements en deux ordres consécutifs et doit être plus petite qu'un paramètre donné $μ$.

$$ \mu \ge \frac{\| q_n(a_{mp}) - q_{n-1}(a_{mp}) \|}{\| q_n(a_{mp}) - q_0 \|} \qquad (21) $$

La MAN permet ainsi de déterminer la totalité d'une branche non linéaire complexe en appliquant pas à pas la procédure décrite précédemment.

# 5 Résultats et validation numériques

Dans cette partie, nous présentons des exemples de plaques avec divers types d'excitation et de conditions aux limites. Certains de ces exemples sont pris de la littérature pour valider notre approche. Le matériau est en aluminium avec un module d'Young E=70.E9 Pa, la masse volumique $ρ=2778$ kg/m³ et le coefficient de Poisson $v=0,3$. Les plaques sont modélisées par des éléments triangulaires DKT de coque qui ont trois noeuds et cinq degrés de liberté par noeud (u, v, w, θx, θy). Pour une raison de symétrie, seulement un quart de la plaque est discrétisé avec 121 noeuds (c-à-d. 726 ddl pour une harmonique). Basé sur des travaux précédents, le paramètre de tolérance $η = 10^{-4}$, l'ordre de troncature n=15 [7] et le nombre d'harmoniques H=3 [18] sont adoptés au début pour un calcul automatique et optimal. Par la suite, nous avons fait une étude approfondie sur la valeur de ces paramètres avec comme objectif l'amélioration de la qualité de la solution et du temps de calcul.

## 5.1 Plaque non amortie

Afin de valider notre approche, nous considérons la réponse non linéaire d'une plaque carrée complètement encastrée Tableau 1. La longueur de la plaque est égale à 240 fois l'épaisseur, la charge appliquée est supposée adimensionnelle d'une amplitude $P_0^d = 0.2$ (d'après la référence [17], $P_0^d = c.P^d / ρh^2ω_1^2$ avec $c = \iint \phi dxdy / \iint \phi^2 dxdy$ où $\phi$ désigne la forme modale normalisée, $P^d$ l'amplitude de la force extérieur uniformément répartie (N/m²), $ω_1$ la pulsation linéaire de la plaque, $ρ$ la masse volumique et h l'épaisseur de la plaque). Le Tableau 1 montre que nos résultats sont très proches de ceux trouvés dans la littérature.
<table><caption>Tableau 1 - Le rapport de pulsation ω/ω₁ de vibration forcée non linéaire d'une plaque carrée complètement encastrée sous une force harmonique uniformément répartie</caption><thead><tr><th colspan="2">ANM [5]</th><th colspan="2">HFEM [20]</th><th colspan="3">Autres méthodes [15], [16] et [17]</th><th colspan="2">La nouvelle approche (ANM)</th></tr><tr><th>W<sub>c</sub>/h</th><th>ω / ω<sub>1</sub></th><th>W<sub>c</sub>/h</th><th>ω / ω<sub>1</sub></th><th>W<sub>c</sub>/h</th><th colspan="2">ω / ω<sub>1</sub></th><th>W<sub>c</sub>/h</th><th>ω / ω<sub>1</sub></th></tr><tr><th></th><th></th><th></th><th></th><th></th><th>FEM+ LIN</th><th>Elliptique</th><th>Perturbation</th><th></th></tr></thead><tbody><tr><td>0,200</td><td>0,216</td><td>0,200</td><td>0,243</td><td>0,2</td><td>0,103</td><td>0,120</td><td>0,123</td><td>0,201</td></tr><tr><td>-0,200</td><td>1,433</td><td>-0,207</td><td>1,428</td><td>-0,2</td><td>1,418</td><td>1,420</td><td>1,420</td><td>-0,206</td></tr><tr><td>0,400</td><td>0,753</td><td>--</td><td>--</td><td>0,4</td><td>0,737</td><td>0,748</td><td>0,748</td><td>0,404</td></tr><tr><td>-0,400</td><td>1,251</td><td>--</td><td>--</td><td>-0,4</td><td>1,243</td><td>1,249</td><td>1,249</td><td>-0,401</td></tr><tr><td>0,600</td><td>0,893</td><td>0,601</td><td>0,897</td><td>0,6</td><td>0,875</td><td>0,895</td><td>0,896</td><td>0,608</td></tr><tr><td>-0,600</td><td>1,209</td><td>-0,590</td><td>1,212</td><td>-0,6</td><td>1,197</td><td>1,212</td><td>1,212</td><td>-0,601</td></tr><tr><td>0,800</td><td>0,991</td><td>--</td><td>--</td><td>0,8</td><td>0,962</td><td>0,994</td><td>0,995</td><td>0,808</td></tr><tr><td>-0,800</td><td>1,215</td><td>--</td><td>--</td><td>-0,8</td><td>1,194</td><td>1,220</td><td>1,221</td><td>-0,801</td></tr><tr><td>1,000</td><td>1,077</td><td>1,001</td><td>1,080</td><td>1</td><td>1,036</td><td>1,082</td><td>1,085</td><td>1,003</td></tr><tr><td>-1,000</td><td>1,246</td><td>-0,995</td><td>1,248</td><td>-1</td><td>1,214</td><td>1,254</td><td>1,256</td><td>-1,003</td></tr></tbody></table>

$$ (L/h=240, P_0^d = 0.2). $$

## 5.2 Plaques amorties

Dans cette section, nous considérons des plaques élastiques, homogènes et isotropes simplement appuyées dont les caractéristiques géométriques et les coefficients d'amortissement sont donnés dans le Tableau 2. La première plaque (notée 1 dans le Tableau 1) est rectangulaire et est soumise à une excitation harmonique transversale uniformément répartie. La seconde est quant à elle carrée et est soumise à une excitation harmonique transversale concentrée au centre.

<table><caption>Tableau 2 - Caractéristiques géométriques et matériau des plaques étudiées.</caption><thead><tr><th>Plaque N°</th><th>Longueur L (m)</th><th>Largeur l (m)</th><th>épaisseur h (m)</th><th>Coefficients d'amortissement</th></tr></thead><tbody><tr><td>1</td><td>0,6</td><td>0,3</td><td>0.001</td><td>α = 0, β</td></tr><tr><td>2</td><td>0,3</td><td>0,3</td><td>0.001</td><td>α = 0.065, β = 0</td></tr></tbody></table>

Dans la Figure 1, nous schématisons le rapport des amplitudes maximales de vibration de la plaque 1 sur son épaisseur (W<sub>c</sub>/h) en fonction du rapport de la pulsation non linéaire sur la pulsation linéaire de la plaque (ω/ω<sub>1</sub>). La force est harmonique avec des amplitudes égales à 5 N/m² et 10N/m², pour des valeurs de la pulsation autour de la première pulsation de résonance. Nous avons considéré les coefficients d'amortissement β=10<sup>-3</sup> et β=10<sup>-4</sup>. Nous remarquons que l'approche actuelle donne les mêmes résultats que ceux calculés par HFEM [24]. Pour le deuxième exemple, plaque 2, nous avons utilisé les fonctions approximatives (22) pour développer les déplacements de la surface moyenne u, v et w comme l'a fait Amabili [2]. Nous avons appliqué la méthode asymptotique numérique pour obtenir la variation du déplacement par rapport à la pulsation. Après cela, nous avons comparé la première composante du vecteur généralisé de déplacement transversal w<sub>11</sub> à celle donnée par les travaux précédents Figure 2a. Nous remarquons que nos résultats coïncident également parfaitement avec les résultats trouvés dans la littérature. À partir de ces résultats nous avons calculé le déplacement vertical au centre de la plaque Figure 2b. Nous avons comparé ce déplacement obtenu d'une part analytiquement et d'autre part par la méthode des éléments finis. Ces deux méthodes utilisent la méthode asymptotique numérique.
$$
\left\{
\begin{aligned}
u(x, y, t) &= \sum_{i=1}^{m} \sum_{j=1}^{n} u_{2i,j}(t) \sin\left(\frac{2\pi i}{a}x\right) \sin\left(\frac{\pi j}{b}y\right) \\
v(x, y, t) &= \sum_{i=1}^{m} \sum_{j=1}^{n} v_{i,2j}(t) \sin\left(\frac{\pi i}{a}x\right) \sin\left(\frac{2\pi j}{b}y\right) \\
w(x, y, t) &= \sum_{i=1}^{m} \sum_{j=1}^{n} w_{i,j}(t) \sin\left(\frac{\pi i}{a}x\right) \sin\left(\frac{\pi j}{b}y\right)
\end{aligned}
\right.
\qquad (22)
$$

<figure><img src="image_2.png" /><figcaption>Figure 1 - Vibration non linéaire forcée d'une plaque rectangulaire amortie pour différentes excitations et coefficients d'amortissement.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>(a) le premier paramètre généralisé</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(b) déplacement du centre de la plaque</figcaption></figure>

<caption>Figure 2 - Vibration non linéaire forcée d'une plaque carrée amortie F=1.74N, α=0.065, β=0.</caption>

Dans Figure 3 et Figure 4, nous avons comparé les réponses de deux plaques: une est complètement encastrée (CE) et l'autre est simplement appuyée (SA). Nous prenons les forces adimensionnelles comme dans la section 5.1. Nous montrons également la réponse libre, de ces plaques, calculée par la méthode asymptotique numérique. Ces figures montrent clairement que la condition aux limites simplement appuyées apportent une réponse "plus souple" que les plaques encastrées.
<figure><img src="image_5.png" /><figcaption>Figure 3 - Vibration non linéaire forcée d'une plaque carrée amortie
$\alpha=2*0.065*\omega_1$ et $\beta=0$</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4 - Vibration non linéaire forcée d'une plaque rectangulaire amortie
$L/l=2$, $\beta=0.0001$ et $\alpha=0$</figcaption></figure>

En plus des résonances principales, l'approche actuelle nous permet d'obtenir les résonances
super harmoniques. Dans la Figure 5, nous présentons l'influence de la variation du nombre
d'harmonique H. Dans ces essais, nous avons choisis les valeurs suivantes pour les paramètres
de la MAN: n=20, η=10⁻⁴ et μ=10⁻⁴ (voir la section 5.3 pour plus de détails). Ainsi, on remarque
que avec quatre harmoniques, une super harmonique apparaît à ω ≈ 0.27ω₁. Tandis que, avec
six harmoniques, nous avons deux super harmoniques. Une est à ω ≈ 0.25ω₁ et la seconde est à
ω ≈ 0.27ω₁. Il est par conséquent plus juste de tenir compte d'un grand nombre d'harmoniques
dans l'étude de problèmes de vibrations non linéaires de plaques amorties. Cependant, le temps
de calcul augmente avec ce nombre d'harmoniques. Dans la prochaine section, nous étudierons en détail l'évolution du coût de calcul avec l'augmentation du nombre d'harmoniques.

<figure><img src="image_7.png" /><figcaption>Figure 5 - Vibration non linéaire forcée d'une plaque rectangulaire amortie F=40N/m², α=0, β=0.001.</figcaption></figure>

## 5.3 Coût de calcul

Dans cette section, nous avons étudié l'influence des paramètres de la MAN (ordre de troncature n, paramètre de tolérance η), du nombre d'harmoniques H et la tolérance de Padé μ sur le temps de calcul. La force d'excitation est supposée égale à 40N/m² et les coefficients d'amortissement α=0 et β=10⁻⁴. Nous avons fait le calcul autour du premier mode propre c.-à-d. de ω = 0 jusqu'à ω = 2ω₁. D'abord, nous prenons H=3, η = 10⁻⁴ et nous changeons l'ordre de troncature n Tableau 3. Avec l'ordre n=20, on a le plus petit nombre de pas et le meilleur résidu final. Après cela, nous avons pris n=20 et nous avons changé la tolérance Tableau 4. on remarque qu'avec η = 10⁻³ on a le plus petit nombre de pas mais cela conduit à une valeur élevée du résidu final. Cependant, cette qualité de solution peut-être améliorée en ajoutant à la fin du calcul une étape de correction. Finalement, pour ce travail, nous avons choisi une valeur η = 10⁻⁴.

<table><caption>Tableau 3 - Le nombre de pas et le résidu final de la représentation polynomiale à ω = 2ω₁ en fonction du nombre de troncature (n). η = 10⁻⁴, H = 3.</caption><tr><td>Ordre n</td><td>19</td><td>20</td><td>21</td></tr><tr><td>Nombre de pas</td><td>70</td><td>67</td><td>75</td></tr><tr><td>résidu final</td><td>0.11</td><td>0.09</td><td>0.11</td></tr></table>
<table><caption>Tableau 4 - Le nombre de pas et le résidu final de la représentation polynomiale à $\omega = 2\omega_l$ en fonction de la tolérance. n=20, H = 3.</caption><tr><td>Tolérance η</td><td>10<sup>-3</sup></td><td>10<sup>-4</sup></td><td>10<sup>-5</sup></td></tr><tr><td>Nombre de pas</td><td>57</td><td>67</td><td>83</td></tr><tr><td>résidu final</td><td>1.045</td><td>0.090</td><td>0.004</td></tr></table>

Nous choisisons maintenant les paramètres suivants: H=3, n=20 et $\eta = 10^{-4}$. Nous utilisons les approximants de Padé avec plusieurs valeurs de tolérance $\eta$ Tableau 5. La représentation polynomiale demande 67 pas pour avoir la courbe entière de $\omega = 0$ jusqu'à $\omega = 2\omega_l$. On remarque qu'avec $\mu = 10^{-3}$ le nombre de pas avec les approximants de Padé est presque la moitié de celui obtenu avec les séries polynomiales. En outre, les approximants de Padé améliorent la qualité de la solution. Le Tableau 5 montre en effet que le résidu final des représentations rationnelles avec $\mu = 10^{-4}$ et $\mu = 10^{-5}$ est vraiment meilleur que celui qu'on obtient avec les approximations polynomiales.

<table><caption>Tableau 5 - Le nombre de pas et le résidu final de la représentation de padé à $\omega = 2\omega_l$ en fonction de la tolérance de Padé. n=20, $\eta = 10^{-4}$ et H = 3</caption><tr><td>Tolérance μ</td><td>10<sup>-3</sup></td><td>10<sup>-4</sup></td><td>10<sup>-5</sup></td></tr><tr><td>Nombre de pas</td><td>36</td><td>50</td><td>55</td></tr><tr><td>Résidu final</td><td>0.1750</td><td>0.0151</td><td>0.0006</td></tr></table>

Finalement, nous avons calculé le temps d'un pas pour différents nombres d'harmoniques. Dans cet essai, nous avons choisi un maillage composé de 435 noeuds (c.-à-d. 2610 ddl pour une harmonique). Un pas de la méthode asymptotique numérique demande une décomposition de matrice, le calcul des seconds membres et enfin le calcul des approximants de Padé. Cette dernière méthode ne nécessite pas un grand temps de calcul. Pour la décomposition de la matrice tangente, le temps croît avec l'augmentation du nombre d'harmoniques. La matrice tangente est en effet constituée de $(2*H-1)^2$ blocs. Dans cette étude la triangulation de matrice est réalisée par une méthode directe classique (méthode de Crout). On définit alors le temps CPU relatif comme étant le rapport entre le temps de décomposition de l'opérateur tangent $\overline{K}_t$ du problème non linéaire (17) et le temps de décomposition de l'opérateur tangent K dans le problème dynamique linéaire qui contient seulement une harmonique:

$$CPU^{\text{relative}} = \frac{CPU(\overline{K}_t)}{CPU(K)} \qquad (23)$$

Le Tableau 6 montre que pour six harmoniques, le temps de la décomposition de la matrice est 127 fois plus grand que le temps de résolution d'un problème linéaire. Ce temps est très important. On peut cependant remarquer que cette décomposition de matrice doit être également faite avec d'autre méthode de résolution (par exemple la méthode de Newton-Raphson). En outre, la méthode de Newton-Raphson nécessite un grand nombre de pas [1] tandis que notre méthode nous permet d'obtenir un grand tronçon de la solution avec une seule décomposition de matrice.

Le calcul du second membre dépend du nombre d'harmoniques et de l'ordre de troncature. Dans la Figure 6, nous avons tracé l'évolution du rapport, temps de calcul du second membre sur le temps de décomposition d'une matrice, en fonction de nombre de troncature n pour chaque nombre d'harmonique. On remarque qu'en augmentant le nombre d'harmoniques, le temps de calcul du second membre devient négligeable par rapport au temps de la
décomposition.

<table><caption>Tableau 6 - Temps de décomposition d'une matrice non linéaire par rapport à ce de la matrice linéaire</caption><thead><tr><th>Nombre d'harmoniques</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th></tr></thead><tbody><tr><th>CPU<sup>relative</sup></th><td>2,84</td><td>12,61</td><td>34,81</td><td>70,27</td><td>127,05</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>Figure 6 - Temps de calcul du second membre par rapport au temps de la décomposition en fonction de l'ordre de troncature</figcaption></figure>

# 6 Conclusion

Dans ce travail, nous avons développé une méthode asymptotique numérique pour résoudre les problèmes de vibrations non linéaires forcées en tenant compte de l'amortissement structural des plaques. Ces plaques sont soumises à des excitations transversales harmoniques. Une grande partie de la solution non linéaire est obtenue en résolvant une séquence de problèmes linéaires qui ont tous la même matrice de rigidité. L'itération de cette méthode, menant à une technique de continuation nous a permis d'obtenir les courbes non linéaires de résonance en n'importe quel intervalle désiré des amplitudes. Les résultats numériques montrent une bonne correspondance avec les résultats issus de la littérature et ceci pour plusieurs valeurs de l'amortissement. En plus des résonances principales, l'approche actuelle a permis d'obtenir les résonances super harmoniques en prenant plusieurs harmoniques dans l'étude. L'augmentation de ce nombre conduit alors à des matrices de rigidité de grande taille. Le temps de calcul pour obtenir une grande partie de la solution croît également. Néanmoins, le temps de calcul total avec la méthode asymptotique numérique est inférieur à celui obtenu avec un algorithme itératif classique [1].

Dans cet article, nous avons aussi montré que même si le temps de calcul du second membre croît il reste négligeable par rapport au temps de décomposition. Pour réduire le temps de calcul total, on pourrait utiliser un solveur linéaire adapté aux problèmes à plusieurs seconds membres [9], comme c'est le cas de la méthode numérique asymptotique. On pourrait également réduire le nombre de pas en changeant la définition du paramètre de chemin [14]. Ceci fait l'objet des travaux actuels.
Références

[1] F. Abdoun, L. Azrar, E.M. Daya, M. Potier-Ferry. Forced harmonic response of viscoelastic structures by an asymptotic numerical method. Computers and Structures 87, 91–100, 2009

[2] M. Amabili. Nonlinear vibrations of rectangular plates with different boundary conditions: theory and experiments. Computers and Structures, 82, 2587-2605, 2004

[3] M. Amabili. Nonlinear vibrations of circular cylindrical panels. Journal of Sound and Vibration 281, 509–535, 2005

[4] M. Amabili. Theory and experiments for large-amplitude vibrations of rectangular plates with geometric imperfections. Journal of Sound and Vibration 291, 539–565, 2006

[5] L. Azrar, E. H. Boutyour and M. Potier-Ferry. Non-linear forced vibrations of plates by an asymptotic numerical method. Journal of Sound and Vibration, 252(4), 657–674, 2002

[6] J-L. Batoz, G. Dhatt. Modélisation des structures par éléments finis. Vol 3: coques. Hermès, Paris 1992.

[7] E. H. Boutyour, H. Zahrouni, M. Potier-Ferry and M. Boudi. Bifurcation points and bifurcated branches by Asymptotic Numerical Method and Padé approximants. International Journal for Numerical Method in Engineering, 60:1987-2012, 2004

[8] E. H. Boutyour, E. M. Daya and M. Potier-Ferry. A harmonic balance method for the non-linear vibration of viscoelastic shell. Comptes Rendus Mécanique, 33, 68–73, 2006

[9] J.-M. Cadou, M. Potier-Ferry. A solver combining reduced basis and convergence acceleration with applications to non-linear elasticity. Submitted to Communications in Numerical Methods in Engineering

[10] B. Cochelin, N. Damil and M. Potier-Ferry. Asymptotic numerical methods and pade approximants for non-linear elastic structures. International Journal for Numerical methods in Engineering, 37, 1187–1213, 1994

[11] E. M. Daya, L. Azrar and M. Potier-Ferry. An amplitude equation for the non-linear vibration of viscoelastically damped sandwich beams. Journal of Sound and Vibration, 271(3), 789-813, 2003

[12] A. Elhage-Hussein, M Potier-Ferry and N Damil. A numerical continuation method based on Padé approximants. Int. J. Solids and Structures, 37, 6981-7001, 2000

[13] M. Ganapathi, B. P. Patel, P. Boisse and O. Polit. Flexural loss factor of sandwich and laminated composite beams using linear and non linear dynamic analysis. Composites: Part B 30, 245-256, 1999

[14] J.-J. Gervais and H. Sadiky. A new steplength control for continuation with the asymptotic numerical method. IMA Journal of Numerical Analysis 22-2, 207-229, 2002

[15] C. S. Hsu. On the application of elliptic functions in non-linear forced oscillations. Quarterly Applied Mathematics 17, 393–407, 1960

[16] C. Mei. Finite element displacement method for large amplitude free flexural vibrations of beams and plates. Computers and Structures 3, 163–174, 1973
[17] C. Mei and K. Decha-Umphai. A finite element method for non-linear forced vibrations of rectangular plates. American Institute of Aeronautics and Astronautics Journal 23, 1104–1110, 1985

[18] F. Pérignon. Vibrations forcées de structures, élastiques; non linéaire. Doctorat thesis, université Aix-Marseille II, France, 2004

[19] P. Ribeiro and M. Petyt. Non-linear vibration of beams with internal resonance by the hierarchical finite element method. Journal of Sound and vibration, 224(4), 591-624, 1999

[20] P. Ribeiro and M. Petyt. Nonlinear vibration of plates by the hierarchical finite element and continuation methods. International Journal of Mechanical Sciences 41, 437–459, 1999

[21] P. Ribeiro and M. Petyt. Non-linear vibration of composite laminated plates by the hierarchical finite element method. Composite Structures 46, 197-208, 1999

[22] P. Ribeiro. Periodic vibration of plates with large displacements. In: Proceedings of the 42nd AIAA/ASME/ASCE/AHS/ASC structures, structural dynamics, and material conference and exhibit, Seattle, WA, (paper A01-25098), 2001

[23] P. Ribeiro and M. Petyt. Experimental detection of modal interaction in the non-linear vibration of a hinged-hinged beam. Journal of Sound and Vibration 277, 943–954, 2004

[24] P. Ribeiro. Nonlinear vibrations of simply-supported plates by the p-version finite element method. Finite Elements in Analysis and Design 41, 911–924, 2005

[25] C. Touzé and M. Amabili. Nonlinear normal modes for damped geometrically nonlinear systems: Application to reduced-order modelling of harmonically forced structures. Journal of Sound and Vibration 298, 958–981, 2006

Annexe

Dans ce travail, on a considéré des forces harmoniques. Alors le vecteur inconnu est aussi
harmonique. Donc, l'équation du mouvement (7) est écrite sous la forme de (11) où les matrices
M, D et B sont symétriques et diagonales par blocs. Pour trois harmoniques par exemple,
nous avons:

$$
\bar{B}_1 = \begin{bmatrix} B_1 & 0 & 0 & 0 & 0 \\ 0 & B_1 & 0 & 0 & 0 \\ 0 & 0 & B_1 & 0 & 0 \\ 0 & 0 & 0 & B_1 & 0 \\ 0 & 0 & 0 & 0 & B_1 \end{bmatrix}, \quad \bar{D} = \begin{bmatrix} D & 0 & 0 & 0 & 0 \\ 0 & D & 0 & 0 & 0 \\ 0 & 0 & D & 0 & 0 \\ 0 & 0 & 0 & D & 0 \\ 0 & 0 & 0 & 0 & D \end{bmatrix}
$$

$$
\bar{M} = \begin{bmatrix} 0 & 0 & 0 & 0 & 0 \\ 0 & M & 0 & 0 & 0 \\ 0 & 0 & M & 0 & 0 \\ 0 & 0 & 0 & 4M & 0 \\ 0 & 0 & 0 & 0 & 4M \end{bmatrix}, \quad \bar{C} = \begin{bmatrix} 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & C & 0 & 0 \\ 0 & -C & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 2C \\ 0 & 0 & 0 & -2C & 0 \end{bmatrix}
$$
La matrice $B_{nl}$ consiste en un produit de deux vecteurs A et B développés sous la forme des séries (10). Par exemple, ce produit $A(t) \times B(t)$ avec $H=3$, est comme suit:

$$
\begin{aligned}
AB ={}& a^0 b^0 + \frac{1}{2} (a_c^1 b_c^1 + a_s^1 b_s^1 + a_c^2 b_c^2 + a_s^2 b_s^2) + \\
      & + \left( a^0 b_c^1 + a_c^1 b^0 + \frac{1}{2} (a_c^1 b_c^2 + a_s^1 b_s^2 + a_c^2 b_c^1 + a_s^2 b_c^1) \right) \cos(\omega t) + \\
      & + \left( a^0 b_s^1 + a_s^1 b^0 + \frac{1}{2} (a_c^1 b_s^2 - a_s^1 b_c^2 - a_c^2 b_s^1 + a_s^2 b_c^1) \right) \sin(\omega t) + \\
      & + \left( a^0 b_c^2 + a_c^2 b^0 + \frac{1}{2} (a_c^1 b_c^1 - a_s^1 b_s^1) \right) \cos(2\omega t) + \\
      & + \left( a^0 b_s^2 + a_s^2 b^0 + \frac{1}{2} (a_c^1 b_s^1 + a_s^1 b_c^1) \right) \sin(2\omega t)
\end{aligned}
$$

Pour employer des vecteurs sous la forme des série (10), nous pouvons réécrire ce produit comme suit:

$$
\begin{bmatrix} (AB)^0 \\ (AB)_c^1 \\ (AB)_s^1 \\ (AB)_c^2 \\ (AB)_s^2 \end{bmatrix} =
\begin{bmatrix}
a^0 & \frac{1}{2}a_c^1 & \frac{1}{2}a_s^1 & \frac{1}{2}a_c^2 & \frac{1}{2}a_s^2 \\
a_c^1 & a^0 + \frac{1}{2}a_c^2 & \frac{1}{2}a_s^2 & \frac{1}{2}a_c^1 & \frac{1}{2}a_s^1 \\
a_s^1 & \frac{1}{2}a_s^2 & a^0 - \frac{1}{2}a_c^2 & -\frac{1}{2}a_s^1 & \frac{1}{2}a_c^1 \\
a_c^2 & \frac{1}{2}a_c^1 & -\frac{1}{2}a_s^1 & a^0 & 0 \\
a_s^2 & \frac{1}{2}a_s^1 & \frac{1}{2}a_c^1 & 0 & a^0
\end{bmatrix}
\times
\begin{bmatrix}
b^0 \\
b_c^1 \\
b_s^1 \\
b_c^2 \\
b_s^2
\end{bmatrix}
$$

L'insertion de cette valeur dans l'équation de mouvement (7) et l'équilibrage des équations selon les harmoniques mènent à réécrire $\bar{B}_{nl}(\bar{U})$ sous la forme:

$$
\bar{B}_{nl}(\bar{U}) =
\begin{bmatrix}
\mathbf{B}_{nl}(U^0) & \mathbf{B}_{nl}\left(\frac{1}{2}U_c^1\right) & \mathbf{B}_{nl}\left(\frac{1}{2}U_s^1\right) & \mathbf{B}_{nl}\left(\frac{1}{2}U_c^2\right) & \mathbf{B}_{nl}\left(\frac{1}{2}U_s^2\right) \\
\mathbf{B}_{nl}(U_c^1) & \mathbf{B}_{nl}\left(U^0 + \frac{1}{2}U_c^2\right) & \mathbf{B}_{nl}\left(\frac{1}{2}U_s^2\right) & \mathbf{B}_{nl}\left(\frac{1}{2}U_c^1\right) & \mathbf{B}_{nl}\left(\frac{1}{2}U_s^1\right) \\
\mathbf{B}_{nl}(U_s^1) & \mathbf{B}_{nl}\left(\frac{1}{2}U_s^2\right) & \mathbf{B}_{nl}\left(U^0 - \frac{1}{2}U_c^2\right) & \mathbf{B}_{nl}\left(-\frac{1}{2}U_s^1\right) & \mathbf{B}_{nl}\left(\frac{1}{2}U_c^1\right) \\
\mathbf{B}_{nl}(U_c^2) & \mathbf{B}_{nl}\left(\frac{1}{2}U_c^1\right) & \mathbf{B}_{nl}\left(-\frac{1}{2}U_s^1\right) & \mathbf{B}_{nl}(U^0) & 0 \\
\mathbf{B}_{nl}(U_s^2) & \mathbf{B}_{nl}\left(\frac{1}{2}U_s^1\right) & \mathbf{B}_{nl}\left(\frac{1}{2}U_c^1\right) & 0 & \mathbf{B}_{nl}(U^0)
\end{bmatrix}
$$

Ici le vecteur inconnu (10) est classé selon des blocs d'harmoniques i (i=0, H-1). Après discrétisation, afin de réduire le nombre d'éléments stockés dans les matrices, numériquement on stocke toutes les harmoniques du même degré de liberté dans le même bloc comme suit:

$$
q = \left[ \underbrace{q_{c1}^0, q_{c1}^1, q_{s1}^1, \dots, q_{c1}^i, q_{s1}^i, \dots, q_{s1}^{H-1}}_{i=1} q_{c2}^0, \dots, q_{s2}^{H-1}, \dots, \underbrace{q_{cj}^0, \dots, q_{sj}^{H-1}}_{j=1} q_{cm}^0, \dots, q_{sm}^{H-1} \right]^T
$$

où j est le degré de liberté (j=1, m).