v2.02.018 SDLL18 - Modes de vibration d’un tuyau droit et mince encastré aux deux extrémités

Résumé:

Ce test consiste à rechercher les fréquences propres et les modes de vibration associés à une tuyauterie droite.

Il permet de valider les modélisations éléments finis DKT (QUAD4), TUYAU_3M (SEG3) et TUYAU_6M (SEG3).

Les résultats obtenus sont comparés à une solution de référence analytique.

Solution de référence

Méthode de calcul utilisée pour la solution de référence

Tout d’abord, il faut remarquer qu’il n’existe pas de formule analytique exacte pour le calcul de modes propres dans le contexte de théorie de coques (Love, Reissner, Flügge, etc.).

Cependant, au cours des années, plusieurs méthodes ont été développées pour approcher les modes propres d’une coque cylindrique ayant plusieurs types de conditions de bord. Parmi elles, on y trouve:

• la méthode GDQ (Generalized differential quadrature) s’inspirant de la

procédure de Ritz par Loy et al. (1997) [1],

• la méthode de propagation d’ondes (ou wave propagation) utilisant de

développements de Fourier par Zhang et al. (2001) [2] et par Xuebin (2008) [3],

• la méthode RRM (reverberation-ray matrix) par Tang et al. (2016) [4].

Résultats de référence

La grande majorité des auteurs présentent leurs résultats en utilisant la fréquence normalisée \(\Omega=\omega R((1-\nu^2)\rho/E)^{1/2}\) où \(\omega\) est la vitesse angulaire en \(\textrm{rad}.\textrm{s}^{-1}\), \(R\) est le rayon moyen du tuyau en \(\textrm{m}\), \(\nu\) est le coefficient de Poisson, \(\rho\) est la masse volumique en \(\textrm{kg}.\textrm{m}^{-3}\) et \(E\) est le module de Young en \(\textrm{Pa}\). La fréquence propre est ainsi calculée par \(f=\omega/(2\pi)\).

Comme l’explique Cammalleri et al. (2016) [5], la déformation en vibration libre d’une coque cylindrique consiste en un nombre \(n\) d’ondes sur la section orthogonale à l’axe, et un nombre \(m\) de demi-ondes sur la section contenant l’axe, appelées respectivement ondes circonférentielles et demi-ondes longitudinales. Ainsi, chaque forme modale est caractérisée par une paire de valeurs \(n\) et \(m\) (voir, par exemple, les Figs. 1, 2 et 3).

La forme des ondes circonférentielles est indépendante des conditions aux limites, alors que la forme des demi-ondes longitudinales dépend des conditions aux limites et est similaire aux vibrations de flexion de poutres soumises aux mêmes contraintes.

Figure 2.2-1: Mode de vibration \(n=2\) et \(m=1\)

Figure 2.2-2: Mode de vibration \(n=3\) et \(m=1\)

Figure 2.2-3: Mode de vibration \(n=4\) et \(m=2\)

Le tableau suivant présente la fréquence normalisée du premier mode axial (\(m=1\)) en fonction du mode circonférentiel (\(n=1,\ldots,10\)), telle que rapportée par différents auteurs.

n	m	Loy, 1997	Zhang, 2001	Xuebin, 2008	Tang, 2016
1	1	0,032885	0,034879	0,034879	0,032792
2	1	0,013932	0,014052	0,014052	0,013902
3	1	0,022672	0,022725	0,022726	0,022667
4	1	0,042208	0,042271	0,042272	0,042208
5	1	0,068046	0,068116	0,068116	0,068046
6	1	0,099748	0,099823	0,099823	0,099749
7	1	0,137249	0,137328	0,137329	0,137251
8	1	0,180535	0,180617	0,180618	0,180536
9	1	0,229599	0,229684	0,229684	0,229684
10	1	0,284439	0,284526	0,284527	0,284526

Tableau 2.2-1: Fréquences normalisées pour différents modes \(n\) avec un mode \(m\) fixe

En utilisant les données du matériau et les résultats de référence précédents, les fréquences propres sont:

n	m	Loy, 1997	Zhang, 2001	Xuebin, 2008	Tang, 2016
1	1	28,377367	30,098044	30,098044	28,297115
2	1	12,022304	12,125856	12,125856	11,996417
3	1	19,56429	19,610025	19,610888	19,559975
4	1	36,422439	36,476804	36,477667	36,422439
5	1	58,718757	58,779162	58,779162	58,718757
6	1	86,075281	86,140001	86,140001	86,076144
7	1	118,435921	118,504093	118,504956	118,437647
8	1	155,788597	155,859357	155,860219	155,789459
9	1	198,127266	198,200615	198,200615	198,200615
10	1	245,450204	245,525279	245,526142	245,525279

Tableau 2.2-2: Fréquences propres pour différents modes \(n\) avec un mode \(m\) fixe

Il est important de noter que, bien qu’en général, pour un \(n\) fixe, la fréquence associée à la paire \((n,m+1)\) soit supérieure à celle du mode \((n,m)\), ce n’est pas nécessairement le cas lorsque l’on garde \(m\) fixe et que l’on fait varier \(n\).

Par exemple, il se peut que la fréquence correspondant au mode \((n=2,m=1)\) soit inférieure à celle du mode \((n=1,m=1)\), ce qui peut conduire à des graphiques de fréquences modales qui semblent inhabituels ou contre-intuitifs (voir par exemple Fig. 2.2-4). Ce phénomène dépend fortement de la géométrie du cylindre et des conditions aux limites, ce qui peut rendre l’interprétation des résultats plus complexe.

Figure 2.2-4: Fréquences propres en fonction du mode circonférentiel pour plusieurs modes axiaux

Comme le montrent les résultats de Xuebin (2008) [3], les modes propres par ordre croissant pour le problème de référence sont les suivants :

Ordre	n	m	Xuebin, 2008
1	2	1	12,13
2	3	1	19,61
3	3	2	23,28
4	2	2	28,06
5	1	1	30,09
6	3	3	31,97
7	4	1	36,48
8	4	2	37,38
9	4	3	39,77

Tableau 2.2-3: Les premiers modes propres en ordre croissant

Références bibliographiques

1. Loy CT, Lam KY, Shu C. Analysis of cylindrical shells using generalized differential quadrature. Shock and vibration. 1997;4(3):193-8.

2. Zhang XM, Liu GR, Lam KY. Frequency analysis of cylindrical panels using a wave propagation approach. Applied Acoustics. 2001 May 1;62(5):527-43.

3. Xuebin L. Study on free vibration analysis of circular cylindrical shells using wave propagation. Journal of sound and vibration. 2008 Apr 8;311(3-5):667-82.

4. Tang D, Wu G, Yao X, Wang C. Free Vibration Analysis of Circular Cylindrical Shells with Arbitrary Boundary Conditions by the Method of Reverberation-Ray Matrix. Shock and Vibration. 2016;2016(1):3814693.

5. Cammalleri M, Costanza A. A closed-form solution for natural frequencies of thin-walled cylinders with clamped edges. International Journal of Mechanical Sciences. 2016 May 1;110:116-26.

Post traitement des résultats du modèle tuyau

import numpy as np

def toroidal_to_cartesien(array_in:np.ndarray):

This function takes an array of points in toroidal coordinates and converts them to Cartesian coordinates. The toroidal coordinates are expected to be in the format [u, v, w, phi], where: - u is the axial coordinate, - v is the orthoradial coordinate, - w is the radial coordinate, - phi is the angle described by v w.r.t. the Cartesien coordinates.

The conversion formulas are: - X = u - Y = w*cos(phi) - v*sin(phi) - Z = w*sin(phi) + v*cos(phi)

Parameters: array_in (np.ndarray): A 2D array where each row represents a point in toroidal coordinates [u, v, w, phi].

Returns: np.ndarray: A 2D array where each row represents a point in Cartesian coordinates [X, Y, Z].

radial = array_in[:, 2] cos_phi = np.cos(array_in[:, 3]) sin_phi = np.sin(array_in[:, 3])

array_out = np.copy(array_in) # array_out[:, 0] = array_in[:, 0] array_out[:, 1] = radial*cos_phi - orthor*sin_phi array_out[:, 2] = radial*sin_phi + orthor*cos_phi return array_out[:, :3] # Return only X Y Z

def postprocess_shell_displacement(U_s:np.ndarray, nb_modes:int, angles:np.ndarray):

The input array U_s is expected to be an array containing the displacement components of the middle surface of a shell structure. The structure of U_s is as follows: - For each Fourier mode (from 1 to nb_modes), the displacements are listed in the order: [u^i_m, v^i_m, w^i_m, u^o_m, v^o_m, w^o_m] - The displacement for mode 0 are listed at the end in the following order: [w^o, w^i_1, w^o_1] Here, i denotes the in-plane Fourier displacement and o denotes the out-plane Fourier displacement, and m denotes the Fourier mode. The number of modes, also called nb_modes, is either 3 or 6.

The angles array contains the azimuthal angles at which the displacements are evaluated. The function returns a 2D array where each row corresponds to a different angle. The first three columns are the axial, ortho-radial, and radial displacements, respectively.

block_list = [] for m in range(2, nb_modes+1):

block = np.array([

[np.cos(m*phi), 0, 0], [0, np.sin(m*phi), 0], [0, 0, np.cos(m*phi)]

]) block_list.append(block)

for m in range(2, nb_modes+1):

block = np.array([

[np.sin(m*phi), 0, 0], [0, np.cos(m*phi), 0], [0, 0, np.sin(m*phi)]

]) block_list.append(block)

block = np.array([

[0, 0, 0], [0, np.sin(phi), -np.cos(phi)], [1, np.cos(phi), np.sin(phi)],

]) block_list.append(block) displacement[i, :-1] = np.block(block_list) @ U_s

return displacement

def postprocess_beam_displacement(U_p:np.ndarray, radius:float, angles:np.ndarray):

The input array U_p is expected to be an array containing the displacement components of the section of a beam structure. The structure of U_p is as follows: - The beam displacements [DX, DY, DZ] are the first three components of U_p, - The beam rotations [RX, RY, RZ] are the last three components of U_p.

The angles array contains the azimuthal angles at which the displacements are evaluated. The radius float is the value of the middle radius of the cylindrical shell. The function returns a 2D array where each row corresponds to a different angle. The first three columns are the axial, ortho-radial, and radial displacements, respectively.

displacement[i, 0] = U_p[0] + U_p[5]*radius*np.sin(phi) - U_p[4]*radius*np.cos(phi) displacement[i, 1] = U_p[2]*np.sin(phi) - U_p[1]*np.cos(phi) - U_p[3]*radius displacement[i, 2] = -(U_p[2]*np.cos(phi) + U_p[1]*np.sin(phi))

return displacement

def sum_toroidal_coords(array_list:list):

This function takes a list of arrays, each representing toroidal coordinates, and returns a single array containing the sum of the input arrays. The input arrays are expected to have the same shape, with the last axis representing the angle values.

Parameters: array_list (list): A list of arrays, where each array represents a set of toroidal coordinates.

Returns: np.ndarray: A single array containing the sum of the input arrays, with the same shape as the individual arrays.

if len(array_list) > 1:

for array in array_list[1:]:

array_out[…, :-1] += array[…, :-1]

return array_out

# Define global variables # For importing and exporting import os SUFIX = “6m” # “3m” or “6m” DIRPATH = os.path.dirname(os.path.abspath(__file__)) COORDS = np.load(f »{DIRPATH}/output/coords_{SUFIX}.npy ») DEPLS = np.load(f »{DIRPATH}/output/displacements_{SUFIX}.npy »)

# For post-processing RADIUS = 1 ANGLES = np.linspace(0, 2*np.pi, 120)

# For plotting OFFSET = 2.5 VMIN, VMAX = 0, 2 NBTICKS = 11 NBCOLORS = 2*NBTICKS-1

# Loop through each frequency from matplotlib import colors, cm, gridspec, pyplot as plt nb_freqs, nb_coords = np.shape(DEPLS)[:2] nb_freqs = np.min([nb_freqs, 18]) for ii in range(nb_freqs):

fig = plt.figure(figsize=(9, 4)) gs = gridspec.GridSpec(1, 2, width_ratios=[1, 2.5], figure=fig) ax1 = fig.add_subplot(gs[0]) ax2 = fig.add_subplot(gs[1], projection=”3d”)

# Loop through each coordinate XX_all, YY_all, ZZ_all, VV_all = [], [], [], [] for jj in range(nb_coords):

# Compute the total displacement beam_disp = postprocess_beam_displacement(U_p=DEPLS[ii, jj, :6],

radius=RADIUS, angles=ANGLES)

shell_disp = postprocess_shell_displacement(U_s=DEPLS[ii, jj, 6:],

nb_modes=int(SUFIX[0]), angles=ANGLES)

total_disp = sum_toroidal_coords([beam_disp, shell_disp])

# Define undeformed cylinder in toroidal coordinates old_cylinder = np.array([[COORDS[jj, 0], 0, RADIUS, phi] for phi in ANGLES])

# Compute deformed cylinder in cartesien coordinates new_cylinder = toroidal_to_cartesien(sum_toroidal_coords([old_cylinder, total_disp]))

# Collect the data for later plotting XX_all.append(new_cylinder[:, 0]) YY_all.append(new_cylinder[:, 1]) ZZ_all.append(new_cylinder[:, 2]) VV_all.append(np.linalg.norm(new_cylinder[:, 1:], axis=1)/RADIUS)

# Convert collected data to arrays for plot XX_all = np.array(XX_all) YY_all = np.array(YY_all) ZZ_all = np.array(ZZ_all) VV_all = np.array(VV_all)

# Create the scatter plot for this frequency ax1.scatter(YY_all, ZZ_all, c=VV_all,

norm=colors.Normalize(vmin=VMIN, vmax=VMAX), cmap=cm.get_cmap(“viridis”, NBCOLORS))

ax1.grid(False) ax1.set_xlabel(“Y”) ax1.set_ylabel(“Z”) ax1.set_xlim([-OFFSET, OFFSET]) ax1.set_ylim([-OFFSET, OFFSET]) ax1.set_aspect(“equal”)

sc = ax2.scatter(XX_all, YY_all, ZZ_all, c=VV_all,

norm=colors.Normalize(vmin=VMIN, vmax=VMAX), cmap=cm.get_cmap(“viridis”, NBCOLORS))

ax2.grid(False) ax2.set_xlabel(“X”) ax2.set_ylabel(“Y”) ax2.set_zlabel(“Z”) ax2.set_xlim([0, 20]) ax2.set_ylim([-OFFSET, OFFSET]) ax2.set_zlim([-OFFSET, OFFSET]) ax2.set_xticks([0, 10, 20])

cbar = fig.colorbar(sc, ax=ax2, shrink=0.8, pad=0.1) cbar.set_ticks(np.linspace(VMIN, VMAX, NBTICKS)) cbar.set_label(“Gonflement normalisé”)

fig.tight_layout() fig.savefig(f »{DIRPATH}/output/tuyau_{SUFIX}_mode_{str(ii+1)}.png ») plt.close(fig)

Modélisation A

Caractéristiques de la modélisation

Figure 3.1-1

Caractéristiques du maillage

Nombre de nœuds: 12663

Nombre de mailles et types: 12600 QUAD4, 326 SEG2

Grandeurs testées et résultats

Incertitudes sur les résultats en prenant les solutions de référence:

	n	m	Ref	DKT	error (%)
1	2	1	12	11,99	0,022
2	3	1	19,56	19,56	0,007
3	3	2	23,28	23,09	0,814
4	2	2	28,06	27,16	3,221
5	1	1	28,3	28,3	0,002
6	3	3	31,97	31,46	1,585
7	4	1	36,42	36,42	0,001
8	4	2	37,38	37,27	0,292

Tableau 3.3-1

Modélisation B

Caractéristiques de la modélisation

Figure 4.1-1

Caractéristiques du maillage

Nombre de nœuds: 961

Nombre de mailles et types: 480 SEG3

Grandeurs testées et résultats

Incertitudes sur les résultats en prenant la solution analytique de référence:

	n	m	Ref	TUYAU 3M	error (%)
1	2	1	11,996	11,81	1,55
2	3	1	19,56	18,7	4,39
3	3	2	23,28	22,32	4,11
4	2	2	28,06	27,04	3,63
5	1	1	28,3	30,2	-6,71
6	3	3	31,97	30,84	3,54
7	4	1	36,42	–	–
8	4	2	37,38	–	–

Tableau 4.3-1

Remarque

Par construction, un élément TUYAU_3M n’utilise que les trois premiers modes de Fourier et n’est donc pas capable de calculer les fréquences de vibration qui correspondent à un mode circonférentiel supérieur à 3. Cela explique pourquoi il n’y a pas de résultats dans les deux dernières lignes du tableau 4.3-1.

Modélisation C

Caractéristiques de la modélisation

Idem que la modélisation B.

Caractéristiques du maillage

Nombre de nœuds: 121

Nombre de mailles et types: 60 SEG3

Grandeurs testées et résultats

Incertitudes sur les résultats en prenant la solution de référence:

	n	m	Ref	TUYAU 6M	error (%)
1	2	1	11,996	11,82	1,46
2	3	1	19,56	18,79	3,92
3	3	2	23,28	22,41	3,73
4	2	2	28,06	27,05	3,59
5	1	1	28,3	30,2	-6,71
6	3	3	31,97	30,92	3,29
7	4	1	36,42	35,03	3,81
8	4	2	37,38	35,88	4,01

Tableau 5.3-1

Synthèse des résultats

Les résultats obtenus avec les modélisations TUYAU_3M (SEG3) et TUYAU_6M (SEG3) sont satisfaisants. L’écart observé est inférieur à 5% pour la plupart des modes.

Les résultats obtenus avec la modélisation DKT (QUAD4) sont aussi satisfaisants. L’écart observé est inférieur à 2% pour la plupart des modes.