u2.05.05 Calcul de structure en fatigue vibratoire#

Contexte industriel#

Ce document a pour but de décrire la mise en œuvre d’un calcul de structure en fatigue vibratoire. On s’intéresse plus particulièrement aux études sur les machines tournantes (ailettes de turbines par exemple). Ces composants sont sollicités par un chargement statique (effort centrifuge lié à la rotation de la machine) et par un chargement dynamique (vibrations induites par l’écoulement du fluide). La durée de vie de la structure dépend à la fois de la partie dynamique (amplitude des variations des contraintes) et de la partie statique (contrainte moyenne).

Le chargement statique est généralement bien connu; la partie dynamique est quant à elle difficile tant à mesurer qu’à calculer. On ne peut donc pas mettre directement en œuvre les outils classiques de calcul de la durée de vie, basés sur l’évolution des contraintes au cours du temps (confer remarque ci-dessous).

A l’inverse, il peut être intéressant d’estimer l’amplitude de variation maximale admissible par la structure, pour un ou plusieurs modes propres de vibration. Cette amplitude peut ensuite être comparée avec les vibrations mesurées sur site (mesures BVM par exemple).

La démarche d’une telle étude est la suivante:

calcul de la contrainte liée au chargement statique ;
calcul de la contrainte associée au mode propre considéré ;
application des critères de fatigue pour calculer, en chaque nœud du maillage, la vibration maximale admissible pour avoir une endurance illimitée de la structure.

La dernière étape peut se faire avec l’opérateur CALC_FATIGUE(TYPE_CALCUL = “FATIGUE_VIBR”). Le principe du calcul, sa mise en œuvre dans Code_Aster et un exemple industriel sont décrits dans les chapitres suivants.

Remarques :

Un certain nombre d’hypothèses sont introduites pour la calcul, cf. paragraphe suivant. Il s’agit ici d’estimer, de manière conservative, un niveau vibratoire acceptable. Si “évolution temporelle des contraintes dans la structure est connue, il faut utiliser les fonctionnalités classiques des opérateurs CALC_FATIGUE ou POST_FATIGUE(avec des critères uniaxiaux ou multiaxiaux, des méthodes de comptage de cycles, …).

En dehors des turbines, la méthodologie présentée ici pourrait s’appliquer également à l’analyse vibratoire de lignes de tuyauterie (la contrainte moyenne correspondant alors à la pression interne), ou pour les éoliennes.

Le calcul nécessite de connaître au préalable la contrainte à la rupture et la limite d’endurance du matériau.

Présentation de la méthode de calcul#

Cadre du calcul et hypothèses#

On se place dans le cadre d’une décomposition de la contrainte dans la structure sur une base modale :

\({\sigma}_{\mathrm{total}}(t)={\sigma}_{\mathrm{stat}}+\sum_{i=1}^{\infty}{\alpha}_{i}{\sigma}_{\mathrm{mod}}^{i}\cos({\omega}_{i}t+{\phi }_{i})\) ,

en notant \({\sigma}_{\mathrm{stat}}\) la contrainte liée au chargement statique et \({\sigma}_{\mathrm{mod}}^{i}\) les contraintes associées aux modes propres de la structure. \({\omega}_{i}\) et \({\phi }_{i}\) sont respectivement la pulsation (connue) et le déphasage (inconnu) du mode i.

Le calcul classique consisterait à évaluer le dommage, les contributions modales et les déphasages étant connus. Ici, réciproquement, on cherche à estimer les contributions modales maximales associées à une endurance illimitée de la structure.

Il n’est pas possible d’identifier de manière générale la contribution maximale de chacun des modes. On est donc conduit à introduire trois hypothèses simplificatrices : critère de fatigue uniaxial ; imposition a priori du poids relatif des modes ; maximisation de l’amplitude.

Critère de fatigue uniaxial : l’utilisation d’un critère de fatigue uniaxial (méthode de Wöhler) revient à supposer que les directions principales du chargement statique et du chargement dynamique sont les mêmes. Cette hypothèse semble licite pour les structures usuelles visées (ailettes, lignes de tuyauterie, …) ; elle induit un conservatisme sans doute excessif dans le cas général. Dans toute la suite, la notation \(\sigma\) correspondra à la norme de la contrainte (von Mises ou Tresca).

A noter qu’une approche multiaxiale (avec des critères de type Crossland ou Dang Van) serait également possible. Ces critères nécessitent toutefois de connaître, en plus de la limite d’endurance en traction du matériau, la limite d’endurance en cisaillement pur. Cette donnée est souvent indisponible pour les matériaux des centrales, et la fonctionnalité aurait été alors difficilement utilisable.

Poids relatif des modes : on suppose que le poids relatif des différents modes propres considérés est connu. Ce poids relatif peut par exemple être estimé à partir de mesures sur site. On ne considère également qu’un nombre limité de modes propres \(N\) (dans la pratique, on aura le plus souvent \(N<4\) ou 5). Autrement dit : \({\sigma}_{\mathrm{total}}(t)={\sigma}_{\mathrm{stat}}+\alpha \sum_{i=1}^{N}{\beta}_{i}{\sigma}_{\mathrm{mod}}^{i}\cos({\omega}_{i}t+{\phi }_{i})\) .

Le coefficient \(\alpha\) est le paramètre que l’on cherche à calculer.

Maximisation de l’amplitude : dans l’expression de la contrainte ci-dessus, les déphasages \({\phi }_{i}\) sont inconnus. La contrainte alternée \({S}_{\mathrm{alt}}\) , définie comme la demi-amplitude de variation de la contrainte , est alors calculée comme suit: \({S}_{\mathrm{alt}}=\alpha \sum_{i=1}^{N}{\beta}_{i}{\sigma}_{\mathrm{mod}}^{i}\) . Cette définition de \({S}_{\mathrm{alt}}\) est conservative.

Calcul des vibrations admissibles#

Sous les hypothèses introduites ci-dessus, le calcul se résume à identifier le coefficient \(\alpha\) associé à une endurance illimitée de la structure.

La contrainte statique \({\sigma}_{\mathrm{stat}}\) correspond ici à une contrainte moyenne, généralement non nulle. Il est nécessaire de prendre en compte cette contrainte dans la courbe de fatigue de Wöhler. Cette prise en compte se fait classiquement à l’aide du diagramme de Haigh [r7.04.01], soit avec la droite de Goodman, soit avec la parabole de Gerber.

En notant \({S}_{l}\) la limite d’endurance et \({S}_{u}\) la limite à la rupture du matériau, on a :

avec la droite de Goodman :

\({S}_{\mathit{alt}}^{\max}={S}_{l}(1-\frac{{\sigma}_{\mathit{stat}}}{{S}_{u}})\) , soit \(\alpha ={S}_{l}(1-\frac{{\sigma}_{\mathrm{stat}}}{{S}_{u}})/\sum_{i=1}^{N}{\beta}_{i}{\sigma}_{\mathrm{mod}}^{i}\) .

avec la parabole de Gerber :

\({S}_{\mathit{alt}}^{\max}={S}_{l}(1-\frac{{\sigma}_{\mathit{stat}}^{2}}{{S}_{u}^{2}})\) , soit \(\alpha ={S}_{l}(1-\frac{{\sigma}_{\mathit{stat}}^{2}}{{{S}_{u}}^{2}})/\sum_{i=1}^{N}{\beta}_{i}{\sigma}_{\mathit{mod}}^{i}\) .

../../../../_images/100004A600001E6300000E5E50450299226B20251.svg — Fig. 45 Diagramme de Haigh#

Ce calcul est fait pour tous les nœuds ou points de Gauss de la structure. Les zones où \(\alpha\) est le plus faible correspondent aux zones qui limitent la durée de vie de la structure.

Pour passer du coefficient \(\alpha\) à l’amplitude de vibration admissible, une opération supplémentaire est à réaliser. Cette opération est décrite dans la u2.05.05-interpretation-resultats.

Mise en œuvre du calcul#

Calcul des contraintes statiques et dynamiques#

L’estimation du niveau vibratoire admissible nécessite le calcul au préalable des contraintes statiques et dynamiques :

le résultat correspondant au chargement statique (effort centrifuge, pression, …) peut être calculé avec MECA_STATIQUE ouSTAT_NON_LINE. Un seul pas de temps doit être présent dans la structure de données résultat.Le calcul des contraintes \({\sigma}_{\mathrm{stat}}\) à partir du champ de déplacement se fait à l’aide de la commande CALC_CHAMP(option SIEQ_ELGA ou SIEQ_ELNO selon que l’on souhaite calculer le dommage aux noeuds ou au points de Gauss).
Le calcul des \(N\) modes propres considérés peut se faire avec la commande CALC_MODES. Le calcul des contraintes modales \({\sigma}_{\mathrm{mod}}^{i}\) se fait à l’aide de la commande CALC_CHAMP(option SIEQ_ELGA ou SIEQ_ELNO selon que l’on souhaite calculer le dommage aux noeuds ou au points de Gauss).

Définition des propriétés matériau#

Deux paramètres matériau sont nécessaires pour le calcul :

la valeur de la limite à la rupture du matériau \({S}_{u}\) . Ce paramètre doit être introduit dans l’opérateur DEFI_MATERIAU [u4.43.01] (mot clé facteur RCCM, opérande Su).
la limite d’endurance du matériau \({S}_{l}\) . Ce paramètre correspond au premier point de la courbe de Wöhler (opérateur DEFI_MATERIAU, mot clé FATIGUE, opérande WOHLER).

Appel à CALC_FATIGUE#

L’utilisation de CALC_FATIGUE (TYPE_CALCUL = “FATIGUE_VIBR”) nécessite de renseigner les paramètres suivants :

les contraintes statiques et modales préalablement calculées ;
les propriétés matériau ;
la liste des modes à considérer (NUME_MODE) et leur poids relatif (FACT_PARTICI), correspondant aux coefficients \({({\beta}_{i})}_{1\le i\le N}\) ;
le choix du mode de prise en compte des contraintes moyennes (Gerber ou Goodman : opérande CORR_SIGM_MOYE) ;
le choix du lieu de calcul du dommage : aux nœuds (OPTION=”DOMA_ELNO_SIGM”) ou aux points de Gauss (OPTION=”DOMA_ELGA_SIGM”).

Interprétation des résultats#

En sortie de l’opérateur CALC_FATIGUE, on dispose d’un champ (aux nœuds ou aux points de Gauss) de la valeur admissible de \(\alpha\) : les zones où \(\alpha\) est le plus faible correspondent aux zones qui limitent la durée de vie de la structure.

Pour passer du coefficient \(\alpha\) à l’amplitude de vibration admissible, une opération supplémentaire est à réaliser. Supposons par exemple qu’on s’intéresse à l’amplitude des déplacement en un point donné, qu’on notera \(\partial \tilde{u}\) . Ce point peut correspondre par exemple à la position d’un capteur, où à la zone d’amplitude de vibration maximale.

On note \({\tilde{u}}_{\mathrm{mod}}^{i}\) le déplacement au point d’intérêt associé au mode i. L’amplitude de vibration admissible au point d’intérêt est alors :

\(\partial \tilde{u}=\min(\alpha )\sum_{i=1}^{N}{\beta}_{i}{\tilde{u}}_{\mathrm{mod}}^{i}\)

Ce calcul est illustré dans l’exemple ci-dessous.

La valeur minimale de \(\alpha\) peut être obtenue de trois manières différentes :

soit en post-traitant dans Aster le champ résultat (POST_RELEVE_T, OPERATION=”EXTREMA”);
soit en visualisant le champ résultat ;
soit dans le fichier message. L’information est imprimée comme suit :

Amplitude de vibration maximale admissible par la structure : 1.318438

Remarques et conseils#

Quelques remarques peuvent être formulées sur cette fonctionnalité :

Le résultat du calcul est nécessairement très sensible à la qualité du calcul des contraintes. L’utilisateur doit donc s’assurer que la finesse de son maillage est suffisante. L’utilisation d’éléments quadratiques est indispensable.
De même, il est conseillé de faire le calcul aux points de Gauss, l’interpolation pour le passage points de Gauss – nœuds étant source d’imprécision.
Les zones critiques correspondent le plus souvent aux singularités géométriques. L’utilisateur devra vérifier si la singularité correspond à une réalité physique ou à une discrétisation trop grossière. A noter que l’utilisation du remaillage (logiciel Homard, opérateur MACR_ADAP_MAIL) ne permet pas de converger vers une solution stable non nulle si la géométrie initiale est discrétisée trop grossièrement.

Exemple#

Description du cas de calcul#

Le post-traitement en fatigue vibratoire va être illustré sur le cas des ailettes des aubages 5-12 des turbines basse pression des centrales REP. Les ailettes sont fixées au disque par des doigts et des broches. Le diamètre extérieur des aubages est de l’ordre de 2 mètres.

Les propriétés matériaux retenues pour cette étude sont les suivantes :

limite d’endurance : \({S}_{l}=400\mathit{MPa}\)
limite à la rupture : \({S}_{u}=800\mathrm{MPa}\) .

Le chargement statique correspond à la rotation de la roue à \(1500\mathit{tours}/min\) . Trois des premiers modes correspondent respectivement à la flexion tangentielle (\(f=100\mathit{Hz}\) ), à la flexion axiale (\(f=\mathrm{180 }\mathit{Hz}\) ) et à la torsion (\(f=360\mathit{Hz}\) )

Remarque: Les propriétés matériau et les chargements utilisés ici ne correspondent pas au cas réel.

Fig. 46 Maillage d’un paquet d’ailette de l’aubage 5-12#

Résultats obtenus#

On s’intéresse ici à l’amplitude maximale de vibration admissible en tête d’ailette (groupe de nœuds “NO_CAPT” , correspondant à la position sur site du capteur BVM), afin que la durée de vie dans le premier doigt de la première ailette du paquet d’aubage soit infinie (groupe de mailles appelé “DOI_1_1”).

Résultats issus de CALC_FATIGUE:

L’appel à CALC_FATIGUE, pour un calcul du dommage aux points de Gauss pour le premier mode de la structure, se fait comme suit :

DMG_MOD1= CALC_FATIGUE( TYPE_CALCUL = “FATIGUE_VIBR”,

OPTION = “DOMA_ELGA_SIGM”,

HISTOIRE = _F( RESULTAT = PRECONT,

MODE_MECA = MODE,

NUME_MODE=1,

FACT_PARTICI = 1,),

DOMMAGE = “WOHLER”,

CORR_SIGM_MOYE = “GOODMAN”,

MATER = MATE,

)

Puis, pour extraire le coefficient \(\alpha\) minimale sur la zone d’intérêt :

IMPR_RESU(FORMAT=”RESULTAT”, RESU=_F(CHAM_GD=DMG_MOD1,GROUP_MA=”DOI_1_1”),);

Les résultats obtenus, pour les trois premiers modes de la structure et pour une combinaison entre les deux premiers modes sont donnés dans le tableau suivant.

	Mode 1 (flexion tangentielle)	Mode 2 (flexion axiale)	Mode 3 (torsion)	Mode 1 + 0,5 Mode 2
Contrainte statique maximale	685 MPa
Contrainte modale maximale	13 MPa	51 MPa	27 MPa
\({\alpha}_{\mathrm{Gerber}}^{\min}\)	1.8	0.8	1.3	1.5
\({\alpha}_{\mathrm{Goodman}}^{\min}\)	3.3	1.5	2.3	2.8

Utilisation des résultats :

Une opération supplémentaire est à réaliser pour passer du coefficient \(\alpha\) à l’amplitude de vibration admissible en un nœud donné. Le nœud considéré ici correspond à la position d’un des capteurs de déplacement.

Cette opération peut par exemple être fait par des boucles python simples :

# Paramètres du calcul (les mêmes que dans CALC_FATIGUE)

fact_partici = [1., 0.5]

nume_mode = [1,2]

nbmode = len(nume_mode)

# EXTRACTION DES DEPLACEMENTS DE LA TABLE

TDEP = POST_RELEVE_T(ACTION=_F( GROUP_NO=”NO_CAPT”,

INTITULE=”Depl. modal”,

RESULTAT=MODE,

NOM_CHAM=”DEPL”,

TOUT_CMP=”OUI”,

OPERATION=”EXTRACTION”,) ,) ;

DX = [None]*nbmode

DY = [None]*nbmode

DZ = [None]*nbmode

tdepl = TDEP.EXTR_TABLE()

for i in range(nbmode) :

dmod = tdepl.NUME_MODE==nume_mode[i]

DX[i] = dmod[“DX”].values()[“DX”]

DY[i] = dmod[“DY”].values()[“DY”]

DZ[i] = dmod[“DZ”].values()[“DZ”]

# EXTRACTION DU COEFFICIENT AMIN

TDMG12=POST_RELEVE_T(ACTION=_F( GROUP_NO=”ELTS”,

INTITULE=”(MODE 1) + 0,5 * (MODE 2)”,

CHAM_GD=D_MOD12,

TOUT_CMP=”OUI”,

OPERATION=”EXTREMA”,) ,) ;

AMIN = TDMG12.EXTR_TABLE()

AMIN = AMIN.EXTREMA==”MIN”

AMIN = AMIN[“VALE”].values()[“VALE”]

# CALCUL DU DEPLACEMENT MAXIMAL ADMISSIBLE

DXmax = 0.

DYmax = 0.

DZmax = 0.

for i in range(nbmode) :

DXmax = DXmax + AMIN[0] * fact_partici[i] * abs(DX[i][0])

DYmax = DYmax + AMIN[0] * fact_partici[i] * abs(DY[i][0])

DZmax = DZmax + AMIN[0] * fact_partici[i] * abs(DZ[i][0])

Dtot = sqrt(DXmax*DXmax + DYmax*DYmax +DZmax*DZmax ).

Cet exemple peut être retrouvé dans le cas test SDLV129A. Les résultats pour la présente étude sont donnés dans le tableau suivant.

Cas	\(\mathit{DX}\) (\(\mathit{mm}\) )	\(\mathit{DY}\) (\(\mathit{mm}\) )	\(\mathit{DZ}\) (\(\mathit{mm}\) )	\({\alpha}_{\mathrm{Gerber}}^{\min}\)	\({\mathit{DX}}^{max}\) (\(\mathit{mm}\) )	\({\mathit{DY}}^{max}\) (\(\mathit{mm}\) )	\({\mathit{DZ}}^{max}\) (\(\mathit{mm}\) )	\({D}^{max}\) (\(\mathit{mm}\) )
Mode 1	-0.02	0.07	0.99	1.8	0.04	0.11	1.78	1.79
Mode 2	0.07	0.98	0.21	0.8	0.06	0.78	0.77	0.8
Mode 3	0.02	-0.18	-0.82	1.3	0.03	0.23	1.07	1.09
Mode 1 + 0,5 Mode 2	0.05	0.56	1.1	1.5	0.08	0.83	1.64	1.84