Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Écoulements compressibles 1D

Écoulement 1D

Les équations de bilan pour un écoulement adiabatique 1D d’un gaz parfait s’écrivent

ρt+ρux=0ρut+ρu2x+px=0t(ρe+12ρu2)+(ρe+12ρu2+p)ux=0\begin{aligned} \frac{\partial\rho}{\partial t}+\frac{\partial\rho u}{\partial x} & = & 0\\ \frac{\partial\rho u}{\partial t}+\frac{\partial\rho u^{2}}{\partial x}+\frac{\partial p}{\partial x} & = & 0\\ \frac{\partial}{\partial t}(\rho e+\frac{1}{2}\rho u^{2})+\frac{\partial(\rho e+\frac{1}{2}\rho u^{2}+p)u}{\partial x} & = & 0\end{aligned}

A ces équations on ajoute l’équation d’état des gaz parfaits:

e=1γ1pρe=\frac{1}{\gamma-1}\frac{p}{\rho}

Si l’écoulement est isentropique (i.e. thermodynamiquement réversible), on remplace l’équation de bilan d’énergie par la relation (S=cste)(S=cste) i.e.:

p/ργ=cstep/\rho^{\gamma}=cste

Vitesse de propagation

Considérons tout d’abord la propagation d’une perturbation de pression dpdp dans un tuyau de section constante contenant un gaz au repos. Cette onde sonore se propage avec une célérité cc. L’apparition d’une perturbation de pression engendre une augmentation dudu de la vitesse du gaz, perturbation qui se propage aussi à la célérité c.c. Considérons un volume de contrôle ABAB qui suit le front de la perturbation CC (figure ci dessous)

propagation d’une perturbation

Figure 1:propagation d’une perturbation

Ce volume de contrôle se déplace à la vitesse cc par rapport à un référentiel fixe et l’écoulement relatif peut être considéré comme stationnaire en temps. Donc dans la section BB après le front (fluide au repos), la vitesse relative vaut uB=cu_{B}=-c , la pression pp et la masse volumique ρ,\rho, et dans la section AA avant le front, la vitesse relative du fluide vaut uA=ducu_{A}=du-c, la pression p+dpp+dp et la masse volumique ρ+dρ\rho+d\rho.

L’équation de bilan de la masse intégrée entre AA et BB s’écrit:

(duc)(ρ+dρ)=cρ(du-c)*(\rho+d\rho)=-c*\rho

d’où au premier ordre l’accroissement de vitesse:

du=cdρρdu=c\frac{d\rho}{\rho}

De la même façon l’équation de bilan de la quantité de de mouvement intégrée entre AA et BB s’écrit:

(ρ+dρ)(duc)2+(p+dp)=ρc2+p(\rho+d\rho)(du-c)^{2}+(p+dp)=\rho c^{2}+p

et en ne conservant que les termes au premier ordre

dp+c2dρ=2ρcdudp+c^{2}d\rho=2\rho c\,du

et en remplaçant dudu, on obtient- l’expression de la célérité de l’onde de pression:

c2=(dpdρ)c^{2}=(\frac{dp}{d\rho})

Cette expression de la célérité des ondes est en fait générale (i.e. ne dépend pas du type de gaz), à condition de noter que l’expression trouvée suppose des petites perturbations (i.e. un écoulement isentropique). Dans ce cas, on a pour un gaz :

c2=(pρ)s=cstec^{2}=\left(\frac{\partial p}{\partial\rho}\right)_{s=cste}

Pour un gaz parfait, on a montré que l’écoulement isentropique vérifiait

pργ=cste\frac{p}{\rho^{\gamma}}=cste

ce qui nous donne l’expression de la célérité des ondes dans un gaz parfait

c2=γpρ=γRTc^{2}=\gamma\frac{p}{\rho}=\gamma RT

qui dépend donc de T\sqrt{T}. Pour de l’air à température ambiante T=320KT=320\,K, R=287R=287, γ=1.4\gamma=1.4 la vitesse du son vaut c=360m/sc=360\,m/s

Remarque: pour un écoulement isentropique, on peut remplacer l’équation de bilan de l’énergie par la relation (valable pour n’importe quel gaz):

c2=dpdρc^{2}=\frac{dp}{d\rho}

Ondes sonores

On considère une perturbation de pression p(x,t)p'(x,t) dans un écoulement d’air stationnaire en moyenne (U0,U_{0}, p0,ρ0p_{0},\rho_{0}). Cette perturbation induit une perturbation de masse volumique ρ\rho' et de vitesse uu'. La pression p0+pp_{0}+p' , la masse volumique ρ0+ρ\rho_{0}+\rho' et la vitesse U0+u\textrm{U}_{0}+u' sont solutions des équations d’Euler. De plus les fluctuations étant faibles, l’écoulement est isentropique et on néglige les fluctuations devant les valeurs moyennes U0,ρ0U_{0},\rho_{0} et p0p_{0} en ne conservant que les termes d’ordre 1.

En se plaçant dans un référentiel en translation uniforme U0U_{0}, la fluctuation pp' est solution de l’équation des ondes:

2pt2γp0ρ02px2=0\frac{\partial^{2}p'}{\partial t^{2}}-\frac{\gamma p_{0}}{\rho_{0}}\frac{\partial^{2}p}{\partial x^{2}}=0

La fluctuation de vitesse uu' est définie par:

ut=1ρ0px\frac{\partial u'}{\partial t}=-\frac{1}{\rho_{0}}\frac{\partial p'}{\partial x}

et la fluctuation de masse volumique ρ\rho':

ρ=ρ0γpp0\rho'=\frac{\rho_{0}}{\gamma}\frac{p'}{p_{0}}

Cette équation des ondes décrit la propagation d’une onde de pression p(x,t)p'(x,t) avec une célérité c0=γp0/ρ0c_{0}=\sqrt{\gamma p_{0}/\rho_{0}} . La solution générale de cette équation est de la forme

p(x,t)=F(xc0t)+G(x+c0t)p'(x,t)=F(x-c_{0}t)+G(x+c_{0}t)

Les fonctions FF et GG sont déterminées par les conditions initiales p(x,0)p'(x,0) et pt(x,0)\frac{\partial p'}{\partial t}(x,0)

Ainsi pour une perturbation sinusoïdale p(x,t)=Beiω(xc0t)p'(x,t)=Be^{i\omega(x-c_{0}t)}, la fluctuation de vitesse uu' vaut u(x,t)=1ρ0c0Beiω(xc0t)u'(x,t)=\frac{1}{\rho_{0}c_{0}}Be^{i\omega(x-c_{0}t)}, et la fluctuation de masse volumique ρ(x,t)=1c02Beiω(xc0t)\rho'(x,t)=\frac{1}{c_{0^{2}}}Be^{i\omega(x-c_{0}t)}.

Ordre de grandeurs

Dans de l’air au repos ρ0=1kg/m3\rho_{0}=1\,kg/m^{3} , p0=105Pap_{0}=10^{5}\,Pa , γ=1.4\gamma=1.4 on a c0=370m/sc_{0}=370\,m/s

Les ondes de pression ont une puissance acoustique de l’ordre de 1 à 100 db. La puissance acoustique en décibel est liée à la puissance moyenne ww' de l’onde par la relation:

db=120+10log10(w)db=120+10\,log_{10}(w')

et la puissance moyenne ww' est donné par:

w=ωc02π0ωc02πp(x,t).u(x,t)dt=B2ρ0c0w'=\frac{\omega c_{0}}{2\pi}\int_{0}^{\frac{\omega c_{0}}{2\pi}}p'(x,t).\overline{u'(x,t)}\,dt=\frac{B^{2}}{\rho_{0}c_{0}}

Cela donne l’ordre de grandeurs des amplitudes BB des ondes de pression: B=104PaB=10^{-4}\,Pa à 1Pa1\,Pa, de l’amplitude des fluctuations de vitesse 0.3106m/s0.3\,10^{-6}\,m/s à 0.3102m/s0.3\,10^{-2}\,m/s et de l’amplitude des fluctuations de masse volumique 0.9109kg/m30.9\,10^{-9}\,kg/m^{3} à 0.9105kg/m30.9\,10^{-5}\,kg/m^{3}. On constate donc que les ondes sonores sont quasiment incompressibles. De plus le déplacement ξ\xi des particules est très faible. En effet l’équation de la trajectoire d’une particule s’écrit:

dξdt=u(ξ,t)avecξ(0)=x\frac{d\xi}{dt}=u'(\xi,t)\,\,avec\,\,\xi(0)=x

En linéarisant l’équation, i.e. en approximant u(ξ,t)=u(x,t)u'(\xi,t)=u'(x,t) , on intègre cette équation:

ξ=iBρ0ωc02eiω(xc0t)=12πfiBρ0c0eiω(xc0t)\xi=\frac{iB}{\rho_{0}\omega c_{0}^{2}}e^{i\omega(x-c_{0}t)}=\frac{1}{2\pi f}\frac{iB}{\rho_{0}c_{0}}e^{i\omega(x-c_{0}t)}

en notant f=ωc02πf=\frac{\omega c_{0}}{2\pi} la fréquence temporelle de l’onde. Pour des ondes acoustiques ( f=10hzf=10\,hz à 10khz10\,khz), avec f=330hzf=330\,hz on obtient une amplitude du déplacement de 0.4109m0.4\,10^{-9}\,m à 0.4105m0.4\,10^{-5}\,m.

Exemple de propagation d’une onde sonore

propagation d’une onde sonore

Figure 2:propagation d’une onde sonore

animation de la propagation d’une onde sonore

Figure 3:animation de la propagation d’une onde sonore

Cône de Mach

On considère l’émission de perturbations (ondes sonores) par une source. Différents cas se présentent suivant que la source se déplace à une vitesse supérieure ou inférieure à la vitesse de propagation des ondes: i.e. la vitesse du son. (voir animations)

  1. Si la source est immobile U0=0U_{0}=0, les surfaces d’ondes sont des sphères

animation de la propagation d’une onde sonore Ma=0

Figure 4:animation de la propagation d’une onde sonore Ma=0Ma=0

  1. Si U0<c0U_{0}<c_{0} (écoulement subsonique)

animation de la propagation d’une onde sonore Ma=0.8

Figure 5:animation de la propagation d’une onde sonore Ma=0.8Ma=0.8

  1. Si U0=c0U_{0}=c_{0} (écoulement transsonique)

animation de la propagation d’une onde sonore Ma=1

Figure 6:animation de la propagation d’une onde sonore Ma=1Ma=1

  1. si U0>c0U_{0}>c_{0} (écoulement supersonique)

animation de la propagation d’une onde sonore Ma=2

Figure 7:animation de la propagation d’une onde sonore Ma=2Ma=2

Dans ce dernier cas, la source se déplace plus vite que l’onde, qui reste contenu dans un cône à l’arrière de la source, que l’on appelle cône de Mach. Le demi angle au sommet α\alpha du cône vérifie:

sinα=c0u0=Ma1sin\,\alpha=\frac{c_{0}}{u_{0}}=Ma^{-1}

Cet angle est inversement proportionnel au nombre de Mach Ma=u0c0Ma=\frac{u_{0}}{c_{0}}

Dans l’animation suivante, on modélise le champ de pression généré par un dipôle harmonique se déplaçant à une vitesse supersonique constante. On visualise le champ de pression généré par le dipôle.

propagation sonore supersonique

Figure 8:propagation sonore supersonique

Équations caractéristiques

Les équations d’Euler forment un système hyperbolique, caractéristique d’un phénomène de propagation. Pour étudier les propriétés de ce système, nous considérerons un écoulement compressible instationnaire 1D isentropique. Le système d’équations s’écrit

ρt+ρux=0ρut+ρu2x+px=0\begin{aligned} \frac{\partial\rho}{\partial t}+\frac{\partial\rho u}{\partial x} & = & 0\\ \frac{\partial\rho u}{\partial t}+\frac{\partial\rho u^{2}}{\partial x}+\frac{\partial p}{\partial x} & = & 0\end{aligned}

avec la relation isentropique

p/ργ=cste=Kp/\rho^{\gamma}=cste=K

Nous allons transformer ce système en éliminant pp et ρ,\rho, pour ne conserver que 2 variables uu et c=γpρc=\sqrt{\frac{\gamma p}{\rho}}.

D’après la relation (23), on déduit la relation entre ρ\rho et cc:

c2=Kγργ1 soit log(ρ)=2γ1log(c)+cstec^{2}=K\gamma\rho^{\gamma-1}\,\,\mbox{ soit }\,log(\rho)=\frac{2}{\gamma-1}log(c)+cste

d’où

dρρ=2γ1dcc\frac{d\rho}{\rho}=\frac{2}{\gamma-1}\frac{dc}{c}

En remplaçant dρd\rho dans l’équation de bilan de la masse

ρt+uρx+ρux=0\frac{\partial\rho}{\partial t}+u\frac{\partial\rho}{\partial x}+\rho\frac{\partial u}{\partial x}=0

on obtient une équation sur uu et cc

2γ1ct+2uγ1cx+cux=0\frac{2}{\gamma-1}\frac{\partial c}{\partial t}+\frac{2u}{\gamma-1}\frac{\partial c}{\partial x}+c\frac{\partial u}{\partial x}=0

De même en remplaçant dans l’équation de quantité de mouvement:

ut+uux+1ρpx=0\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+\frac{1}{\rho}\frac{\partial p}{\partial x}=0

dp=c2dρdp=c^{2}d\rho on obtient

ut+uux+2cγ1cx=0\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+\frac{2c}{\gamma-1}\frac{\partial c}{\partial x}=0

En additionnant (27) et (29) et en soustrayant (27) et (29), on obtient les 2 équations caractéristiques:

{t+(u+c)x}(+u+2γ1c)=0{t+(uc)x}(u+2γ1c)=0\begin{aligned} \left\{ \frac{\partial}{\partial t}+(u+c)\frac{\partial}{\partial x}\right\} (+u+\frac{2}{\gamma-1}c) & = & 0\\ \left\{ \frac{\partial}{\partial t}+(u-c)\frac{\partial}{\partial x}\right\} (-u+\frac{2}{\gamma-1}c) & = & 0\end{aligned}

Posons r=u+2cγ1r=u+\frac{2c}{\gamma-1} et s=u+2cγ1s=-u+\frac{2c}{\gamma-1} (invariants de Riemann), le système s’écrit:

rt+(u+c)rx=0st+(uc)sx=0\begin{aligned} \frac{\partial r}{\partial t}+(u+c)\frac{\partial r}{\partial x} & = & 0\\ \frac{\partial s}{\partial t}+(u-c)\frac{\partial s}{\partial x} & = & 0\end{aligned}

dont l’interprétation est la suivante: soit Γ+\Gamma^{+}la famille de courbes de pente dxdt=u+c\frac{dx}{dt}=u+c et Γ\Gamma^{-}celle de pente dxdt=uc\frac{dx}{dt}=u-c, les solutions rr et ss vérifient

  1. r=cster=cste le long des courbes Γ+\Gamma^{+}

  2. s=cstes=cste le long des courbes Γ\Gamma^{-}

En effet soit ξ\xi l’abscisse curviligne le long de la courbe Γ+\Gamma^{+} de pente dxdt=u+c\frac{dx}{dt}=u+c, on a

drdξ=rtdtdξ+rxdxdξ=(rt+rx(u+c))dtdξ=0\frac{dr}{d\xi}=\frac{\partial r}{\partial t}\frac{dt}{d\xi}+\frac{\partial r}{\partial x}\frac{dx}{d\xi}=\left(\frac{\partial r}{\partial t}+\frac{\partial r}{\partial x}(u+c)\right)\frac{dt}{d\xi}=0

Connaissant les courbes caractéristiques, on peut alors en déduire la solution en tout point.

courbes caractéristiques

Figure 9:courbes caractéristiques

Soit M=(x,t)M=(x,t) un point quelconque du plan (x,t)(x,t), il se trouve à l’intersection de 2 caractéristiques: Γ+\Gamma^{+} issue de A=(xA,0)A=(x_{A},0) et Γ\Gamma^{-} issue de B=(xB,0)B=(x_{B},0). On a donc:

r(M)=r(A) et s(M)=s(B)r(M)=r(A)\,\mbox{ et }\,s(M)=s(B)

Connaissant rr et ss on en déduit l’état du fluide en MM.

Cette démarche attractive possède cependant un inconvénient majeur: on ne sait pas en général déterminer les courbes caractéristiques Γ+\Gamma^{+} et Γ\Gamma^{-}, car elles dépendent de la solution ! Cependant cette approche nous informe sur les propriétés des équations d’Euler pour un écoulement de gaz:

  1. ces équations traduisent une propagation (équations hyperboliques)

  2. les vitesses de propagation sont u+cu+c et ucu-c ( et uu dans le cas général)

  3. la vitesse de propagation n’est pas forcement constante et dépend de la solution (problème non linéaire)

  4. pour les petites perturbations (ondes acoustiques) ucu\ll c et la vitesse de propagation vaut ±c\pm c

  5. le sens de propagation change suivant le signe de ucu-c

  6. certaines quantités (invariants de Riemann) se conservent le long des caractéristiques

  7. ces équations peuvent générées des discontinuités (chocs), comme on le verra dans la suite

Création des discontinuités

Considérons un écoulement d’air se déplaçant à la vitesse 2U02U_{0} dans un milieu initialement au repos, correspondant au champ de vitesse initial continu:

u(x,0)=U0(1tanhxl)u(x,0)=U_{0}(1-\tanh\frac{x}{l})

La répartition de pression et de masse volumique à l’instant initial est telle que la célérité du son c(x,t)c(x,t) vérifie:

c(x,0)=γ12u(x,0)+c0c(x,0)=\frac{\gamma-1}{2}u(x,0)+c_{0}

i.e. on a pour x>2lx>2l un fluide au repos, avec une pression p0p_{0} et une masse volumique ρ0\rho_{0} t.q. c02=γp0ρ0c_{0}^{2}=\frac{\gamma\,p_{0}}{\rho_{0}}, pour x<2lx<-2l un fluide en mouvement à la vitesse 2U02U_{0} (avec 2U0<c02U_{0}<c_{0} écoulement initial subsonique), avec une pression p1p_{1} et une masse volumique ρ1\rho_{1} t.q. (c0+(γ1)U0)2=γp1ρ1(c_{0}+(\gamma-1)U_{0})^{2}=\frac{\gamma\,p_{1}}{\rho_{1}}, et entre 2l<x<2l-2l<x<2l une transition continue en tangente hyperbolique tanhxl\tanh\frac{x}{l} (rq tanh20.96\tanh2\simeq0.96 ).

Avec ces conditions initiales, les invariants de Riemann ont pour valeurs:

s(x,0)=2γ1c(x,0)u(x,0)=2γ1c0, et r(x,0)=2U0(1tanhxl)+2γ1c0s(x,0)=\frac{2}{\gamma-1}c(x,0)-u(x,0)=\frac{2}{\gamma-1}c_{0},\,\mbox{\, et\, }\,r(x,0)=2U_{0}(1-\tanh\frac{x}{l})+\frac{2}{\gamma-1}c_{0}
schéma des courbes caractéristiques

Figure 10:schéma des courbes caractéristiques

D’après ce qui précède, ss est constant le long des courbes caractéristiques Γ\Gamma^{-} qui sont issue de l’axe OxOx . Or ss est constant à l’instant initial (i.e. sur l’axe OxOx ), donc ss doit être constant dans tout le plan (x,t)(x,t)

Pour un point M du plan, on a (voir schéma):

s(M)=2γ1c0, et r(M)=2U0(1tanhx0l)+2γ1c0s(M)=\frac{2}{\gamma-1}c_{0},\,\mbox{\, et\, }\,r(M)=2U_{0}(1-\tanh\frac{x_{0}}{l})+\frac{2}{\gamma-1}c_{0}

Or d’après la définition de rr et ss on a: u=12(rs)u=\frac{1}{2}(r-s) et c=γ12(r+s)c=\frac{\gamma-1}{2}(r+s), ss étant constant dans la plan et rr étant constant le long des courbes Γ+\Gamma^{+}, on en déduit que uu et cc sont forcément constants le long des courbes Γ+\Gamma^{+}. Or les courbes Γ+\Gamma^{+} ont pour pente u+cu+c , ce sont donc des droites. La courbe Γ+\Gamma^{+} issue de AA est la droite d’équation x=(u(xA,0)+c(xA,0))t+xAx=(u(x_{A},0)+c(x_{A},0))*t+x_{A}. La famille des droites Γ+\Gamma^{+} a pour pente:

dtdx=(γ+12U0(1tanhx0l)+c0)1=m01\frac{dt}{dx}=\left(\frac{\gamma+1}{2}U_{0}(1-\tanh\frac{x_{0}}{l})+c_{0}\right)^{-1}=m_{0}^{-1}

Pour x<2lx<-2l , la pente devient constante (((γ+1)U0+c0)1\sim\left((\gamma+1)U_{0}+c_{0}\right)^{-1}) et les caractéristiques sont parallèles, et de même pour x>2lx>2l (la pente vaut c01\sim c_{0}^{-1}). Entre les deux la pente varie continûment et les courbes caractéristiques Γ+\Gamma^{+} doivent se rejoindre au bout d’un temps tst_{s} (intersection des droites issues de x=2lx=-2l et x=+2lx=+2l

ts4l(γ+1)U0t_{s}\simeq\frac{4l}{(\gamma+1)U_{0}}
intersection des courbes caractéristiques

Figure 11:intersection des courbes caractéristiques

Or lorsque des caractéristiques de la même famille se coupent, la solution ne peut plus rester continue, et une discontinuité doit apparaître. L’hypothèse d’écoulement isentropique n’est plus vérifiée et on a apparition d’un choc de compression.

Les courbes caractéristiques ont pour équation (avec les notations (38):

xm0t=x0x-m_{0}t=x_{0}

On a montré que u(x,t)u(x,t) était constant le long de ces caractéristiques, donc u(x,t)=F(xm0t)u(x,t)=F(x-m_{0}t)

Connaissant la solution initiale u(x,t=0)u(x,t=0) (34), on en déduit la solution u(x,t)u(x,t)

u(x,t)=U0(1tanhxm0tl)u(x,t)=U_{0}(1-\tanh\frac{x-m_{0}t}{l})

et la solution c(x,t)c(x,t)

c(x,t)=γ12U0(1tanhxm0tl)+c0c(x,t)=\frac{\gamma-1}{2}U_{0}(1-\tanh\frac{x-m_{0}t}{l})+c_{0}

On note que m0m_{0} (38) vérifie:

m0=u(x,0)+c(x,0)=u(x,t)+c(x,t)m_{0}=u(x,0)+c(x,0)=u(x,t)+c(x,t)

La solution est donc une onde qui se propage à la vitesse absolue m0m_{0} , c’est à dire à la vitesse c(x,t)c(x,t) relativement au fluide. La vitesse de propagation des perturbation par rapport au fluide est donc bien c(x,t)c(x,t), la célérité du son locale. Pour une perturbation faible U0c0U_{0}\ll c_{0}, on retrouve l’approximation de l’acoustique avec des ondes se propageant avec une célérité constante c0c_{0}.

Cette analyse est confirmée par une résolution numérique des équations d’Euler, dont la solution au cours du temps est tracée ci dessous.

ondes de compression

Figure 12:ondes de compression

Si on considère une onde de détente, i.e. une perturbation correspondant à un champ de vitesse opposé au précédent:

u(x,0)=U0(1+tanhxl)u(x,0)=U_{0}(1+\tanh\frac{x}{l})

la perturbation s’évanouit au cours du temps, puisque les courbes caractéristiques divergent (voir figure ci dessous)

ondes de détente

Figure 13:ondes de détente

Écoulement quasi 1D en conduite

On considère un écoulement quasi-1D dans une conduite de section S(x)S(x) variable. On suppose que la variation de section SS est suffisamment faible pour que l’écoulement soit uni-dimensionnel avec une vitesse U(x,t)U(x,t).

écoulement en tuyère

Figure 14:écoulement en tuyère

On considère le volume de contrôle VV de la figure ci-dessus, limité par les sections S(x)S(x) et S(x+dx)S(x+dx)

Les équations de bilan d’un gaz parfait compressible (voir chapitre 1) intégrées sur le volume VV de frontière Γ\Gamma s’écrivent (après utilisation du théorème de la divergence):

V(ρt+div(ρU))dV=VρtdV+ΓρU.ndS=0V(ρUt+div(ρUU+pId))dV=VρUtdV+Γ((ρUU).n+pn)dS=0VEtt+div(UHt)=VEttdV+SHtU.ndS=0\begin{aligned} \int_{V}(\frac{\partial\rho}{\partial t}+div(\rho\overrightarrow{U}))\,dV= & \int_{V}\frac{\partial\rho}{\partial t}dV+\int_{\Gamma}\rho\overrightarrow{U}.\overrightarrow{n}\,dS & =0\\ \int_{V}(\frac{\partial\rho\overrightarrow{U}}{\partial t}+div(\rho\overrightarrow{U}\otimes\overrightarrow{U}+p\,\overline{\overline{Id}}))\,dV= & \int_{V}\frac{\partial\rho\overrightarrow{U}}{\partial t}dV+\int_{\Gamma}((\rho\overrightarrow{U}\otimes\overrightarrow{U}).\overrightarrow{n}+p\overrightarrow{n})\,dS & =0\\ \int_{V}\frac{\partial E_{t}}{\partial t}+div(\overrightarrow{U}H_{t})= & \int_{V}\frac{\partial E_{t}}{\partial t}dV+\int_{S}H_{t}\overrightarrow{U}.\overrightarrow{n}\,dS & =0\end{aligned}

On a noté Et=ρe+12ρU2E_{t}=\rho e+\frac{1}{2}\rho U^{2} l’énergie totale par unité de volume, et Ht=Et+pH_{t}=E_{t}+p l’entalpie totale.

A ces équations, on rajoute l’équation d’état des gaz parfaits:

e=1γ1pρe=\frac{1}{\gamma-1}\frac{p}{\rho}

Le volume VV considéré est égal au premier ordre en x à VSdxV\simeq S\,dx, donc les intégrales de volume d’une fonction f(x)f(x) s’écrivent:

Vf(x)dVf(x)Sdx\int_{V}f(x)\,dV\simeq f(x)\,S\,dx

Les intégrales de surface d’une fonction F(x)F(x) se décomposent en intégrale sur les sections débitantes S(x)S(x) et S(x+dx)S(x+dx) , et sur la frontière pariétale SwS_{w}:

ΓF(x)dS=S(x)F(x)dS+S(x+dx)F(x+dx)dS+SwF(x)dS=F(x)S(x)+F(x+dx)S(x+dx)+SwF(x)dS\begin{aligned} \int_{\Gamma}F(x)\,dS= & \int_{S(x)}F(x)\,dS+\int_{S(x+dx)}F(x+dx)\,dS+\int_{S_{w}}F(x)\,dS\\ = & F(x)S(x)+F(x+dx)S(x+dx)+\int_{S_{w}}F(x)\,dS\end{aligned}

Pour l’équation de bilan de la masse, F(x)=ρU.nF(x)=\rho\overrightarrow{U}.\overrightarrow{n} est nulle sur les frontières pariétales et vaut +ρU(x+dx)+\rho U(x+dx) en S(x+dx)S(x+dx) et ρU(x)-\rho U(x) en S(x)S(x) (la normale n\overrightarrow{n} est toujours dirigée vers l’extérieur du volume de contrôle). On obtient alors:

ρtS(x)dx+ρ(x+dx)U(x+dx)S(x+dx)ρ(x)U(x)S(x)=0\frac{\partial\rho}{\partial t}S(x)dx+\rho(x+dx)U(x+dx)S(x+dx)-\rho(x)U(x)S(x)=0

soit au premier ordre en x:

Sρt+x(ρUS)=0S\frac{\partial\rho}{\partial t}+\frac{\partial}{\partial x}(\rho US)=0

Pour l’équation de bilan de quantité de mouvement, F(x)F(x)

F(x)=((ρUU+pId).n).exF(x)=\left((\rho\overrightarrow{U}\otimes\overrightarrow{U}+p\,\overline{\overline{Id}}).\overrightarrow{n}\right).\overrightarrow{e_{x}}

est la composante suivant OxOx qui est non nulle sur les frontières pariétales et vaut pnxp\,n_{x}nx=n.exn_{x}=\overrightarrow{n}.\overrightarrow{e_{x}} est la composante suivant OxOx de la normale à la frontière. L’intégrale sur la paroi s’écrit au premier ordre en x:

SwpnxdSp(x+dx)+p(x)2SwnxdS=p(x+dx)+p(x)2(S(x)S(x+dx))\int_{S_{w}}p\,n_{x}dS\simeq\frac{p(x+dx)+p(x)}{2}\int_{S_{w}}n_{x}dS=\frac{p(x+dx)+p(x)}{2}(S(x)-S(x+dx))

Sur S(x)S(x), F(x)F(x) vaut ρU2p-\rho U^{2}-p et sur S(x+dx)S(x+dx) F(x+dx)F(x+dx) vaut +ρU2+p+\rho U^{2}+p, et l’équation de bilan de la quantité de mouvement s’écrit:

ρUtS(x)dx+(ρU2+p)S(x+dx)(ρU2+p)S(x)+p(x+dx)+p(x)2(S(x)S(x+dx))=0\frac{\partial\rho U}{\partial t}S(x)\,dx+(\rho U^{2}+p)S(x+dx)-(\rho U^{2}+p)S(x)+\frac{p(x+dx)+p(x)}{2}(S(x)-S(x+dx))=0

ce qui donne au premier ordre en x:

ρUtS(x)+x(ρU2S)+x(pS)pSx=0\frac{\partial\rho U}{\partial t}S(x)+\frac{\partial}{\partial x}(\rho U^{2}S)+\frac{\partial}{\partial x}(pS)-p\frac{\partial S}{\partial x}=0

soit:

SρUt+x(ρU2S)+Spx=0S\frac{\partial\rho U}{\partial t}+\frac{\partial}{\partial x}(\rho U^{2}S)+S\frac{\partial p}{\partial x}=0

Enfin l’équation de bilan de l’énergie, F(x)=HtU.nF(x)=H_{t}\overrightarrow{U}.\overrightarrow{n} qui s’annulle sur les parois et vaut HtU-H_{t}U en S(x)S(x) et +HtU+H_{t}U en S(x+dx)S(x+dx), ce qui donne au final:

EttSdx+(HtU)S(x+dx)(HtU)S(x)=0\frac{\partial E_{t}}{\partial t}S\,dx+(H_{t}U)S(x+dx)-(H_{t}U)S(x)=0

soit à l’ordre 1 en x:

SEtt+x(HtUS)=0S\frac{\partial E_{t}}{\partial t}+\frac{\partial}{\partial x}(H_{t}U\,S)=0

Les équations de bilan pour un gaz parfait dans une conduite de section S(x)S(x) s’écrivent donc:

Sρt+x(ρUS)=0SρUt+x(ρU2S)+Spx=0SEtt+x(HtUS)=0\begin{aligned} S\frac{\partial\rho}{\partial t}+\frac{\partial}{\partial x}(\rho US) & = & 0\\ S\frac{\partial\rho U}{\partial t}+\frac{\partial}{\partial x}(\rho U^{2}S)+S\frac{\partial p}{\partial x} & = & 0\\ S\frac{\partial E_{t}}{\partial t}+\frac{\partial}{\partial x}(H_{t}U\,S) & = & 0\end{aligned}

avec l’équation d’état: Ht=Et+p=γγ1p+12ρU2H_{t}=E_{t}+p=\frac{\gamma}{\gamma-1}\,p+\frac{1}{2}\rho U^{2}

Écoulement stationnaire en conduite

Dans le cas d’un écoulement stationnaire, le système d’équations (50)- (57) se simplifie.

L’équation de bilan de la masse (50) s’écrit:

ρUS=cste\rho US=cste

i.e. le débit massique (i.e. en kg/s) se conserve .

L’équation de bilan de l’énergie (57) s’écrit:

HtUS=csteH_{t}US=cste

En notant ht=Htρh_{t}=\frac{H_{t}}{\rho} l’entalpie totale par unité de masse, on obtient:

ht=csteh_{t}=cste

i.e. l’entalpie totale par unité de masse se conserve.

Pour l’équation de bilan de la quantité de mouvement (55), on obtient (en utilisant (50):

x(ρU2S)+Spx=ρUSUx+Ux(ρUS)+Spx0=ρUSUx+Spx\begin{aligned} \frac{\partial}{\partial x}(\rho U^{2}S)+S\frac{\partial p}{\partial x} & = & \rho US\frac{\partial U}{\partial x}+U\frac{\partial}{\partial x}(\rho US)+S\frac{\partial p}{\partial x}\\ 0 & = & \rho US\frac{\partial U}{\partial x}+S\frac{\partial p}{\partial x}\end{aligned}

ce qui se simplifie

ρUUx+px=0\rho U\frac{\partial U}{\partial x}+\frac{\partial p}{\partial x}=0

Sous forme différentielle, on peut l’écrire: ρUdU+dp=0\rho U\,dU+dp=0.

Le système d’équation pour un écoulement permanent en conduite s’écrit

ρUS=csteht=csteρUdU+dp=0\begin{aligned} \rho US & = & cste\\ h_{t} & = & cste\\ \rho U\,dU+dp & = & 0\end{aligned}

auquel on adjoint l’équation d’état ht=γγ1pρ+12U2h_{t}=\frac{\gamma}{\gamma-1}\frac{p}{\rho}+\frac{1}{2}U^{2}

Dans le cas où l’écoulement ne présente pas de discontinuité (choc), l’écoulement est aussi isentropique (exercice).

En prenant le logarithme de l’équation (64), on obtient:

logρ+logU+logS=cste\log\rho+\log U+\log S=cste

soit en différenciant:

dρρ+dUU+dSS=0\frac{d\rho}{\rho}+\frac{dU}{U}+\frac{dS}{S}=0

En utilisant l’équation (64) pour remplacer dUdU, on a:

dSS=dρρ+dpρU2=dpρU2(1U2(dpdρ))\frac{dS}{S}=-\frac{d\rho}{\rho}+\frac{dp}{\rho U^{2}}=\frac{dp}{\rho U^{2}}\left(1-\frac{U^{2}}{(\frac{dp}{d\rho})}\right)

En introduisant la célérité du son c=dpdρc=\sqrt{\frac{dp}{d\rho}}, et le nombre de Mach local M=UcM=\frac{U}{c}, l’équation précédente s’écrit:

dSS=dpρU2(1M2)\frac{dS}{S}=\frac{dp}{\rho U^{2}}(1-M^{2})

Interprétons cette relation et ces conséquences si on veut accélérer un gaz:

On note ici une différence importante entre un écoulement subsonique et supersonique:

un convergent accélère un écoulement subsonique, mais décélère un écoulement supersonique.

Tuyère de Laval

Nous allons maintenant étudier le moyen de générer un écoulement supersonique à l’aide d’une tuyère convergente-divergente, dite tuyère de Laval. Ce dispositif est très utilisé pour étudier expérimentalement les écoulements supersoniques.

écoulement en tuyère

Figure 15:écoulement en tuyère

On considère un réservoir, contenant un gaz à la pression p0p_{0} et à la température T0T_{0}, relié à une sortie qui est à une pression pep_{e} et une température TeT_{e} par une conduite convergente divergente de section S(x)S(x) , dont la valeur minimum vaut SS^{*} au col.

Le réservoir est suffisamment grand pour considérer que la vitesse du fluide y est négligeable (U0=0U_{0}=0). Le fluide est quasiment au repos et donc dans des conditions d’arrêt: pression d’arrêt p0p_{0} et température d’arrêt T0T_{0}.

En supposant que l’écoulement est stationnaire et isentropique dans la tuyère, les lois de bilan s’écrivent:

  1. bilan d’énergie (ou premier principe de la thermodynamique):
    l’enthalpie totale ht=h+12U2h_{t}=h+\frac{1}{2}U^{2} se conserve et est égale à l’enthalpie d’arrêt h0h_{0}

    h0=h1+12U12=h2+12U22h_{0}=h_{1}+\frac{1}{2}U_{1}^{2}=h_{2}+\frac{1}{2}U_{2}^{2}
  2. conservation de l’entropie (ou second principe de la thermodynamique):

    s0=s1=s2s_{0}=s_{1}=s_{2}
  3. conservation de la masse:
    le débit massique m˙\dot{m} se conserve dans la tuyère

    ρ1U1S1=ρ2U2S2=m˙=cste\rho_{1}U_{1}S_{1}=\rho_{2}U_{2}S_{2}=\dot{m}=cste
  4. bilan de quantité de mouvement (principe fondamental loi de Newton)

    (p1+ρ1U12)S1+R=(p2+ρ2U22)S2(p_{1}+\rho_{1}U_{1}^{2})S_{1}+R=(p_{2}+\rho_{2}U_{2}^{2})S_{2}

    RR est la résultante des forces de pression sur la paroi entre les 2 sections:

    R=12p.nxdSR=\int_{1}^{2}p.n_{x}\,dS

relations d’arrêt

Nous allons tout d’abord déterminer les relations entre l’état générateur (conditions d’arrêt) et l’état dans une section S1S_{1} quelconque de la conduite (en supposant un écoulement isentropique entre les deux).

La conservation de l’entalpie totale (69) permet d’obtenir une relation entre les températures, en notant que l’enthalpie pour un gaz parfait s’écrit h=CpTh=C_{p}T

h0=CpT0=CpT1+12U12h_{0}=C_{p}T_{0}=C_{p}T_{1}+\frac{1}{2}U_{1}^{2}

soit

T0T1=1+U122CpT1=1+(γ1)U122c12=1+γ12M12\frac{T_{0}}{T_{1}}=1+\frac{U_{1}^{2}}{2C_{p}T_{1}}=1+\frac{(\gamma-1)U_{1}^{2}}{2c_{1}^{2}}=1+\frac{\gamma-1}{2}M_{1}^{2}

en introduisant la vitesse du son c12=γRT1c_{1}^{2}=\gamma RT_{1} , la relation R=CpCvR=C_{p}-C_{v} et le nombre de Mach local M1M_{1}.

En utilisant la relation isentropique pργ=cste\frac{p}{\rho^{\gamma}}=cste , qui s’écrit en fonction de TT: Tργ1=cste\frac{T}{\rho^{\gamma-1}}=cste, on obtiens l’évolution de la masse volumique:

ρ1ρ0=(T1T0)1γ1=(11+γ12M12)1γ1\frac{\rho_{1}}{\rho_{0}}=\left(\frac{T_{1}}{T_{0}}\right)^{\frac{1}{\gamma-1}}=\left(\frac{1}{1+\frac{\gamma-1}{2}M_{1}^{2}}\right)^{\frac{1}{\gamma-1}}

et l’évolution de la pression:

p1p0=(ρ1ρ0)γ=(11+γ12M12)γγ1\frac{p_{1}}{p_{0}}=\left(\frac{\rho_{1}}{\rho_{0}}\right)^{\gamma}=\left(\frac{1}{1+\frac{\gamma-1}{2}M_{1}^{2}}\right)^{\frac{\gamma}{\gamma-1}}

L’évolution de ces quantités est donnée sur la figure ci-dessous:

relation d’arrêt

Figure 16:relation d’arrêt

L’évolution de l’écoulement isentropique à partir de condition d’arrêt correspond donc à une diminution de la température, i.e. à la transformation de l’enthalpie (énergie interne + pression) en énergie cinétique. En effectuant un développement limité en nombre de Mach au voisinage de M1=0M_{1}=0, on obtient les expressions suivantes:

T1T0=1γ12M12p1p0=1γ2M12ρ1ρ0=112M12\begin{aligned} \frac{T_{1}}{T_{0}} & = & 1-\frac{\gamma-1}{2}M_{1}^{2}\\ \frac{p_{1}}{p_{0}} & = & 1-\frac{\gamma}{2}M_{1}^{2}\\ \frac{\rho_{1}}{\rho_{0}} & = & 1-\frac{1}{2}M_{1}^{2}\end{aligned}

On constate que les variations relatives à faible nombre de Mach sont proportionnelles à M12M_{1}^{2}, et sont négligeables pour M1<0.1M_{1}<0.1 (moins de 1%1\%). A température ambiante c=360m/sc=360\,m/s, un nombre de Mach M1=0.1M_{1}=0.1 correspond donc à une vitesse de 36m/s36\,m/s soit 130km/h130\,km/h. Cela justifie l’hypothèse d’écoulement incompressible ρ=cste\rho=cste pour des vitesses telles que M<0.1M<0.1.

ATTENTION: pour un écoulement incompressible, la masse volumique, la pression et la température sont constantes dans l’équation de bilan de masse et d’énergie (les fluctuations sont négligeables). Par contre dans l’équation de bilan de quantité de mouvement le terme en gradient de pression gradp-\overrightarrow{grad}\,p doit être conservé, car il est du même ordre que le terme d’inertie U.gradU\overrightarrow{U}.\overrightarrow{grad}\,\overrightarrow{U}. Ce sont en effet les fluctuations de pression (dites fluctuations dynamiques) qui génèrent l’écoulement. Par contre ces fluctuations sont négligeables devant la pression totale (pression thermodynamique), ce qui permet de considérer les propriétés thermodynamiques d’un fluide incompressible comme constantes.

évolution du Mach dans une tuyère adaptée

Étudions l’évolution de la quantité de mouvement Q=ρUQ=\rho U dans la tuyère. Pour cela exprimons QQ en fonction du nombre de Mach MM , de pp et de TT

Q=pRTMγRT=MpγRTQ=\frac{p}{RT}M\sqrt{\gamma RT}=Mp\sqrt{\frac{\gamma}{RT}}

en remplaçant pp et TT en fonction des conditions d’arrêt (75)- (77)

Q=Mp0γRT0(11+γ12M2)γγ1(1+γ12M2)12Q=Mp_{0}\sqrt{\frac{\gamma}{RT_{0}}}\left(\frac{1}{1+\frac{\gamma-1}{2}M^{2}}\right)^{\frac{\gamma}{\gamma-1}}\left(1+\frac{\gamma-1}{2}M^{2}\right)^{\frac{1}{2}}

soit

Q=p0γRT0M(1+γ12M2)γ+12(γ1)Q=p_{0}\sqrt{\frac{\gamma}{RT_{0}}}\frac{M}{\left(1+\frac{\gamma-1}{2}M^{2}\right)^{\frac{\gamma+1}{2(\gamma-1)}}}

Si la tuyère est adaptée, l’écoulement est sonique au col SS^{*}, et la quantité de mouvement au col QQ^{*} est obtenue en faisant M=1M=1 dans l’expression précédente

Q=p0γRT0(2γ+1)γ+12(γ1)Q^{*}=p_{0}\sqrt{\frac{\gamma}{RT_{0}}}\left(\frac{2}{\gamma+1}\right)^{\frac{\gamma+1}{2(\gamma-1)}}

d’où le rapport QQ=SS\frac{Q^{*}}{Q}=\frac{S}{S^{*}} qui est aussi le rapport des sections d’après la conservation de la masse:

QQ=SS=1M[2γ+1(1+γ12M2)]γ+12(γ1)\frac{Q^{*}}{Q}=\frac{S}{S^{*}}=\frac{1}{M}\left[\frac{2}{\gamma+1}\left(1+\frac{\gamma-1}{2}M^{2}\right)\right]^{\frac{\gamma+1}{2(\gamma-1)}}

Cette relation permet de définir l’évolution du nombre de Mach en fonction de la section de la tuyère.

évolution du nombre de Mach dans une tuyère

Figure 17:évolution du nombre de Mach dans une tuyère

On constate que pour une section donnée, on a deux nombres de Mach possibles, un Mach supersonique ou un Mach subsonique. Donc pour une tuyère de Laval adaptée, on a 2 écoulements possibles. En effet dans la partie convergente, l’écoulement est forcément subsonique (puisqu’on part d’un écoulement au repos). Par contre dans la partie divergente, l’écoulement peut être soit subsonique, soit supersonique suivant la valeur de la pression en sortie de tuyère. En effet ayant la valeur du nombre de Mach, la relation (77) nous fournit la pression correspondante. Sur la figure ci-dessous on a tracé l’évolution des quantités M,p,T,ρM,p,T,\rho pour une tuyère parabolique et pour les 2 types d’écoulements possibles. D’après la relation de bilan de la quantité de mouvement (64), on déduit l’evolution de la vitesse en fonction de celle de la pression, puisque la vitesse augmente (du>0du>0) lorsque la pression diminue dans la tuyère (dp>0dp>0).

évolution du Mach dans une tuyère

Figure 18:évolution du Mach dans une tuyère

évolution de pression, masse volumique et température dans une tuyère

Figure 19:évolution de pression, masse volumique et température dans une tuyère

On constate que si la pression de sortie (en S/S=1.29S/S^{*}=1.29 ,i.e. x=1.7x=1.7) vaut Pe1=0.83P0Pe_{1}=0.83\,P_{0} l’écoulement reste subsonique dans la tuyère et sonique au col, par contre si la pression vaut Pe2=0.22P0Pe_{2}=0.22\,P_{0}, alors l’écoulement est supersonique dans la partie divergente.

Si la pression en sortie PePe est supérieure à Pe1Pe_{1}, alors l’écoulement n’est plus sonique au col et reste subsonique dans la tuyère.

Si la pression en sortie PePe est comprise entre Pe1Pe_{1} et Pe2Pe_{2} , alors l’écoulement n’est plus isentropique et on montrera que dans ce cas un choc se produit dans la tuyère.

Enfin si la pression en sortie est inférieure à Pe2Pe_{2} , alors en sortie l’écoulement n’est plus isentropique et un choc en sortie apparaît.