Bonjour à tous,
Je suis actuellement en Master 2 de mécanique, et je travaille sur un projet où il m’est demandé de reproduire les figures du papier de Tichy (1991) :
Hydrodynamic lubrication theory for the Bingham plastic flow model
(Si vous n’avez pas accès à l’article, je peux vous le transmettre par message privé.)
Je me concentre sur le squeeze film damper, et mon objectif est de retrouver les figures 4, 5 et 6 du papier.
Cependant, malgré de nombreux essais, je n’obtiens pas les résultats attendus.
En particulier, pour les valeurs données par Tichy : τ* = 5 et ε = 0.75, le papier indique q* = -0.3835, tandis que mon code converge vers q* ≈ -0.650.
L’ensemble de mon implémentation est disponible sur GitHub :
http://"https://github.com/KitsuHokk...ricationSolver
Méthode numérique employée
Mon schéma essaie de suivre aussi fidèlement que possible la démarche décrite par Tichy :
- Le problème étant cyclique, je pose le problème en coordonnées angulaires, avec la variable θ = x/R.
- Je crée une grille régulière sur [0, 2π].
Pour chaque point de la grille :
- Je calcule l’épaisseur de film : h = c(1 + ε cos θ).
- Je détermine le gradient de pression newtonien à partir de l’équation (14), afin de connaître le signe du cisaillement (donc l’orientation de τ₀).
- Je construis la cubique de Tichy (équation (26)), dérivée de l’équation (17) dans le cas particulier du squeeze film damper (U₁ = U₂).
- Je calcule les racines de cette équation cubique.
Parmi ces racines :
- je conserve uniquement les racines réelles ;
- je garde celles ayant le même signe que le flux newtonien ;
s’il reste plusieurs racines valides, je sélectionne celle qui est la plus proche de la valeur précédente ou la plus proche du flux newtonien (j’ai testé les deux) afin d’assurer la continuité du profil dp/dx.- J’en déduis les bornes du noyau flottant :
hₐ = (±τ₀ − τ₁)/(dp/dx), h_b = hₐ − 2(±τ₀)/(dp/dx)
où τ₁ est calculé selon l’équation (16a).- Je détermine ensuite le régime local d’écoulement à partir de la condition (24) :
Je sélectionne alors le bon dp/dx selon le régime détecté :
- si 0 < hₐ < h_b < h → floating core ;
- sinon → no core.
- pour no core, j’utilise la formulation newtonienne (équation (14)) ;
- pour floating core, j’utilise la racine filtrée de l’équation (26).
- Une fois tous les dp/dx(θ) obtenus sur la grille, je calcule la pression :
p(θ) = R ∫₀^θ (dp/dx′) dθ′
via une intégration trapézoïdale (fonction cumulative_trapezoid de scipy).- Enfin, j’applique la condition de périodicité :
p(2π; q) − p(0; q) = 0
qui définit le résidu pour la méthode de tir.- Le paramètre q est ensuite mis à jour par une méthode de Newton–Raphson, jusqu’à ce que la périodicité soit respectée.
Observations actuelles
- Le code converge numériquement, mais vers une valeur q différente de celle du papier.
Les profils dp/dx(θ) et p(θ) obtenus sont continus et ressemblent qualitativement aux figures du papier, mais :Lorsque j’utilise la valeur q* = -0.3835 du papier :
- les amplitudes et pentes diffèrent significativement des valeurs publiées ;
- le saut de pente dp/dx observé au passage no core → floating core reste visible, même avec un maillage fin.
- le profil de p(θ) n’est plus périodique ;
- le profil de dp/dx devient continu, mais strictement positif pour chaque θᵢ.
Je termine en précisant que j’ai encore des doutes sur la sélection de la racine dans l’équation (26) et sur le choix du signe de τ₀, au vu des résultats obtenus.
Merci à ceux qui prendront le temps de lire et éventuellement de me donner leur avis.
Cordialement,
David
-----