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.

Couche limite

Marc BUFFAT, Dpt Mécanique, Lyon 1

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rcParams
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 14

Couche limite

from IPython.display import YouTubeVideo

YouTubeVideo('GgVCTNCwfQk', width=600, height=400)
Loading...

Fundamentals of Boundary Layers (MIT)

from IPython.display import IFrame
#IFrame("https://techtv.mit.edu/videos/a081da3d343242bd845ab379f72c281c/embed/",width="560",height="315")

Equation de Blasius

Avec les hypothèses de couche limite bi-dimensionnelle sans gradient de pression (voir Blasius boundary layer), les composantes de vitesse s’écrivent en fonction de ψ(x,y)\psi(x,y), la fonction de courant

u(x,y)=∂ψ∂y    v(x,y)=−∂ψ∂xu(x,y) = \frac{\partial \psi}{\partial y} \;\; v(x,y) = -\frac{\partial \psi}{\partial x}
(1)

dont on cherche une forme auto-similaire en éffectuant un changement de variable

Recherche d’une solution auto-similaire

η=yδ(x)=yU0νx\eta = \frac{y}{\delta(x)} = y \sqrt{\frac{U_0}{\nu x }}
(2)

solution en fonction de courant ψ(x,y)\psi(x,y)

ψ(x,y)=νU0xf(η)\psi(x,y) = \sqrt{\nu U_0 x} f(\eta)
(3)

où f(η)f(\eta) est solution de

2f′′′+ff′′=02 f''' + f f'' = 0
(4)

avec les conditions:

f(0)=0f(0)=0 , f′(0)=0f'(0)=0 et f′′(∞)=1f''(\infty) = 1

Problème: ce n’est pas une ODE classique, car on a pas 3 conditions initiales !!!

Il faut donc trouvez la condition initiale manquante β\beta

f(0)=0 et f′(0)=0 et f′′(0)=βf(0)=0 \mbox{ et }f'(0)=0 \mbox{ et } f''(0)=\beta
(5)

telle que f′′(∞)=1f''(\infty) = 1

→\rightarrow résolution numérique par transformation en ODE d’ordre 1

import sympy as sp
x,y = sp.symbols('x y',positive=True)
nu, U0 = sp.symbols('nu U0',positive=True)
u = sp.Function('u')(x,y)
v = sp.Function('v')(x,y)
p = sp.Function('p')(x)
eta = sp.symbols('eta',positive=True)
f = sp.Function('f')(eta)
psi = sp.Function('psi')(x,y)
eq = u*sp.diff(u,x) + v*sp.diff(u,y) + sp.diff(p,x) - nu*(sp.diff(u,x,x)+sp.diff(u,y,y))
display(eq)
deq = sp.diff(eq,y)
deq1=deq.subs({u:sp.diff(psi,y),v:-sp.diff(psi,x)})
display(deq1)
deq2=deq1.subs(psi,sp.sqrt(nu*U0*x)*f.subs(eta,y*sp.sqrt(U0/x)))
#deq2.simplify().subs(
Loading...
Loading...
import sympy as sp
x,y = sp.symbols('x y')
U0, nu = sp.symbols('U_0 nu')
eta = y*sp.sqrt(U0/(nu*x))
delta = y/eta
f = sp.Function('f')(eta)
psi = sp.sqrt(nu*U0*x)*f
display("psi=",psi)
'psi='
Loading...
U = sp.diff(psi,y)
V = -sp.diff(psi,x)
display(U,V)
Loading...
Loading...

ODE d’ordre 3 transformée en 3 ODEs d’ordre 1

dYdt=F(Y,t)\frac{dY}{dt} = F(Y,t)
(6)

forme canonique: ODE d’ordre 1

dY1dt=Y2,dY2dt=Y3,dY3dt=−12 Y1 Y3\frac{dY_1}{dt} = Y_2 , \frac{dY_2}{dt} = Y_3 , \frac{dY_3}{dt} = -\frac{1}{2}\,Y_1\,Y_3
(7)

avec les C.L.

Y1(0)=0,Y2(0)=0,Y3(0)=β0, t.q. Y∞=Y2(10)=1Y_1(0)=0 , Y_2(0)=0, Y_3(0)=\beta_0, \mbox{ t.q. }Y_\infty = Y_2(10)=1
(8)

### on en déduit la solution autosimilaire fonction de f(η)f(\eta)

η=yU0νx\eta = y \sqrt{\frac{U_0}{\nu x }}
(9)
u(x,y)=U0f′(η)v(x,y)=12νU0x(ηf′(η)−f(η))\begin{align*} u(x,y) &= U_0 f'(\eta) \\ v(x,y) &= \frac{1}{2} \sqrt{\frac{\nu U_0}{x }}(\eta f'(\eta) - f(\eta)) \end{align*}
(10)

résolution numérique

Utilisation de solveur classique ODE

from scipy.integrate import odeint

def ODEBlasius(x,Y):
    '''systeme ODE Blasius'''
    dY = np.array([Y[1],Y[2],-0.5*Y[0]*Y[2]])
    return dY
def Blasius(beta0):
    """resoud BLASIUS pour une valeur de beta"""
    Tmax = 10.
    X  = np.linspace(0,Tmax,500)
    Y0 = np.array([0,0,beta0])
    Y  = odeint(ODEBlasius,Y0,X,tfirst=True)
    return X,Y
def ErrBlasius(beta):
    """calcul l'erreur sur la CL en 0"""
    X,Y = Blasius(beta[0])
    err = Y[-1,1] - 1.0
    #print("beta={} err={}".format(Beta[0],err))
    return err

Problème choix de la bonne valeur de β\beta!

# verification pour beta = 0.3
beta0 = 0.3
X,Y = Blasius(beta0)
print("err=",ErrBlasius([beta0]))
plt.figure(figsize=(10,6))
plt.plot(X,Y[:,0],X,Y[:,1],X,Y[:,2])
plt.legend(["f","f'","f''"])
plt.hlines(1.,0,10.,linestyles='--')
plt.ylim([0,2])
plt.xlabel("$\eta$")
plt.title("solution EDO pour $\\beta$ fixé -> CL # 1");
err= -0.06544374655251373
<Figure size 1000x600 with 1 Axes>

détermination de la solution par optimisation (racine)

from scipy.optimize import fsolve

beta_opt= fsolve(ErrBlasius,[beta0],xtol=1.e-12)[0]
print("beta opt=",beta_opt)
beta opt= 0.3320573347194112

On retrouve bien l’approximation classique de f(η)f(\eta) pour η\eta petit:

f(η)=12βη2 avec β≈0.33205733f(\eta) = \frac{1}{2} \beta \eta^2 \mbox{ avec } \beta \approx 0.33205733
(11)

tracé de la solution

X,Y = Blasius(beta_opt)
plt.figure(figsize=(10,6))
plt.plot(X,Y[:,0],X,Y[:,1],X,Y[:,2])
plt.legend(["f","f'","f''"])
plt.xlabel("$\eta$")
plt.hlines(1.,0,10.,linestyles='--')
plt.ylim(0,1.25)
plt.title("solution autosimilaire");
<Figure size 1000x600 with 1 Axes>
plt.figure(figsize=(8,6))
plt.plot(Y[:,1],X)
plt.ylim(0,5.)
plt.xlim(0,1.)
plt.xlabel("$U/U_0$")
plt.ylabel("$\eta$")
plt.title("Profil de vitesse auto-similaire");
<Figure size 800x600 with 1 Axes>
nu = 15e-6
U0 = 1.0
Xp  = np.array([0.2, 1., 2., 5.,10.])
plt.figure(figsize=(15,6))
for k,x in enumerate(Xp):
    eta0 = np.sqrt(U0/(nu*x))
    plt.subplot(1,len(Xp),k+1)
    ax = plt.gca()
    if k>0: 
        ax.get_yaxis().set_visible(False)
    else:
        plt.ylabel("y")
    plt.ylim(0,0.04)
    plt.xlim(0,1.)
    plt.plot(Y[:,1],(X/eta0))
    plt.title("x={:.2f} ".format(x))
<Figure size 1500x600 with 5 Axes>
# valeur de eta0 tq V=0.99 U0
i0 = np.abs(Y[:,1]-0.99).argmin()
eta0 = X[i0]
print("eta0=",eta0,Y[i0,1])
plt.plot(X,Y[:,1])
plt.hlines(Y[i0,1],X[0],X[-1],linestyles='dashed')
plt.plot([eta0],[Y[i0,1]],'o')
eta0= 4.909819639278557 0.9899968804588336
<Figure size 640x480 with 1 Axes>
Delta = eta0*np.sqrt(nu*Xp/U0)
plt.figure(figsize=(15,6))
plt.plot(Xp,Delta,'-o')
plt.xlabel("x")
plt.ylabel("$\\delta$")
plt.title("Epaisseur de CL : $\\delta =\eta_0\\sqrt{\\nu X/U_0}$");
<Figure size 1500x600 with 1 Axes>

Equation de Falkner Scan

prise en compte d’un gradient de pression avec un champ de vitesse extérieur

ue(x)=U0(xL)mu_e(x) = U_0 \left(\frac{x}{L}\right)^m
(12)

recherche d’une solution auto-similaire

η=yU0(m+1)2νL(xL)(m−1)/2\eta = y \sqrt{\frac{U_0 (m+1)}{2 \nu L}}\left(\frac{x}{L}\right)^{(m-1)/2}
(13)

on se ramène a résoudre une EDO d’ordre 3

f′′′+ff′′+β(1−(f′)2)=0f''' + f f'' + \beta ( 1 - (f')^2) = 0
(14)

avec

β=2mm+1\beta = \frac{ 2m }{m+1}
(15)

où m représente l’angle de la plaque / à l’horizontal

  • m=0m = 0 blasius avec une échelle δ=2νxU0\delta = \sqrt{\frac{2\nu x}{U_0}}

  • m<0m < 0 adverse pressure gradient

  • m>0m > 0 favorable pressure gradient

résolution numérique

def ODEFalkner(x,Y):
    global BETA
    '''systeme ODE Falkner scan'''
    dY = np.array([Y[1],Y[2],-Y[0]*Y[2] + BETA*(Y[1]**2 - 1.0 )])
    return dY
def Falkner(beta0):
    """resoud Falkner pour une valeur de beta"""
    Tmax = 10.
    X  = np.linspace(0,Tmax,50)
    Y0 = np.array([0,0,beta0])
    Y  = odeint(ODEFalkner,Y0,X,tfirst=True)
    return X,Y
def ErrFalkner(beta):
    """calcul l'erreur sur la CL en 0"""
    X,Y = Falkner(beta[0])
    err = Y[-1,1] - 1.0
    #print("beta={} err={}".format(Beta[0],err))
    return err

gradient de pression favorable m > 0

m = 1
BETA = 2*m/(m+1)
beta0 = 1.23
X1,Y1 = Falkner(beta0)
print("err=",ErrFalkner([beta0]))
plt.figure(figsize=(10,6))
plt.plot(X1,Y1[:,0],X1,Y1[:,1],X1,Y1[:,2]);
plt.hlines(1.,0,10.,linestyles='--')
plt.ylim([0,2])
plt.xlabel("$\eta$")
plt.title("solution EDO pour $\\beta$ fixé -> CL # 1");
err= -0.23391875896336956
<Figure size 1000x600 with 1 Axes>

beta_opt= fsolve(ErrFalkner,[beta0])[0]
print("beta opt=",beta_opt)
beta opt= 1.232587684584402

X1,Y1 = Falkner(beta_opt)
plt.figure(figsize=(10,6))
plt.plot(X1,Y1[:,0],X1,Y1[:,1],X1,Y1[:,2])
plt.legend(["f","f'","f''"])
plt.xlabel("$\eta$")
plt.hlines(1.,0,10.,linestyles='--')
plt.ylim(0,2)
plt.title("solution autosimilaire");
<Figure size 1000x600 with 1 Axes>
plt.figure(figsize=(10,6))
plt.plot(Y[:,1],X,'--',label="blasius")
plt.plot(Y1[:,1],X1,label="grad p <0")
plt.ylim(0,5.)
plt.xlim(0,1.)
plt.xlabel("$U/U_0$")
plt.ylabel("$\eta$")
plt.legend()
plt.title("Profil de vitesse auto-similaire");
<Figure size 1000x600 with 1 Axes>

gradient de pression adverse

m = -0.08
BETA = 2*m/(m+1)
beta0 = 0.1
X2,Y2 = Falkner(beta0)
print("err=",ErrFalkner([beta0]))
plt.figure(figsize=(10,6))
plt.plot(X2,Y2[:,0],X2,Y2[:,1],X2,Y2[:,2]);
plt.hlines(1.,0,10.,linestyles='--')
plt.ylim([0,2])
plt.xlabel("$\eta$")
plt.title("solution EDO pour $\\beta$ fixé -> CL # 1");
err= -0.014668533761202918
<Figure size 1000x600 with 1 Axes>
beta_opt= fsolve(ErrFalkner,[beta0])[0]
print("beta opt=",beta_opt)
beta opt= 0.14973564750325108

X2,Y2 = Falkner(beta_opt)
plt.figure(figsize=(10,6))
plt.plot(X2,Y2[:,0],X2,Y2[:,1],X2,Y2[:,2])
plt.legend(["f","f'","f''"])
plt.xlabel("$\eta$")
plt.hlines(1.,0,10.,linestyles='--')
plt.ylim(0,2)
plt.title("solution autosimilaire");
<Figure size 1000x600 with 1 Axes>
plt.plot(Y[:,1],X)
plt.ylim(0,4)
plt.title("Profil de vitesse");
<Figure size 640x480 with 1 Axes>
plt.figure(figsize=(10,6))
plt.plot(Y[:,1],X,'--',label="blasius")
plt.plot(Y1[:,1],X1,label="grad p <0")
plt.plot(Y2[:,1],X2,lw=3,label="grad p >0")
plt.ylim(0,5.)
plt.xlim(0,1.)
plt.xlabel("$U/U_0$")
plt.ylabel("$\eta$")
plt.legend()
plt.title("Profil de vitesse auto-similaire");
<Figure size 1000x600 with 1 Axes>

Analyse dimensionnelle

  1. paramètres : U0U_0 , xx, δ\delta , μ\mu ρ\rho

  2. deux nombres sans dimension: Rex=ρU0xμ≫1Re_x = \frac{\rho U_0 x}{\mu} \gg 1 et ϵ=δx≪1\epsilon = \frac{\delta }{x} \ll 1$

  • bilan de masse

∂u∂x+∂v∂y=0\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0
(16)

ordre de grandeur V0V_0 de v

V0=δxU0=ϵU0V_0 = \frac{\delta}{x} U_0 = \epsilon U_0
(17)

  • bilan de quantité de mouvement suivant x

ρu(∂u∂x+v∂u∂y)=−∂p∂x+μ(∂2u∂x2+∂2u∂y2)\rho u \left(\frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y}\right) = - \frac{\partial p}{\partial x} + \mu \left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right)
(18)

couhe limite laminaire stationnaire sans gradient de pression p≈cstep\approx cste

ρu∂u∂x≈ρv∂u∂y≈μ∂2u∂y2\rho u\frac{\partial u}{\partial x} \approx \rho v\frac{\partial u}{\partial y} \approx \mu \frac{\partial^2 u}{\partial y^2}
(19)

d’où

ρU02x≈μU0δ2\rho \frac{U_0^2}{x} \approx \mu \frac{U_0}{\delta ^2}
(20)

soit

δ2=μxρU0→ϵ≈(Rex)−1/2\delta^2 = \frac{\mu x}{\rho U_0} \rightarrow \epsilon \approx (Re_x) ^{-1/2}
(21)

FIN

Ecoulements de base incompressibles
Approximation de couche limite
Ecoulements de base incompressibles
Approximation de Boussinesq