Résonateur de Helmholtz — Vérification, validation et longueur effective¶
Solveur de différences finies axisymétrique — étude autonome¶
Ce notebook est entièrement autonome : il ne dépend d'aucun fichier local, seulement de
numpy, scipy, matplotlib. Le solveur et toutes les analyses sont écrits dans les
cellules ci-dessous ; tout est recalculé à l'exécution (~3–5 min).
Démarche de vérification & validation (V&V), en trois niveaux distincts :
- vérification du code — le schéma est-il correctement programmé ? → solution manufacturée ;
- vérification de la solution — l'erreur de maillage est-elle bornée ? → convergence GCI ;
- validation — le modèle reproduit-il la réalité ? → mesures publiées.
Sommaire : 1. Théorie · 2. Solveur FDM · 3. Vérification du code (MMS) · 4. Convergence GCI · 5. Validation (Selamet) · 6. Loi d'échelle · 7. Correction de bout (rayonnement) · 8. Pertes & Q · 9. Conclusion.
In [1]:
import numpy as np
import scipy.sparse as sp
from scipy.sparse.linalg import spsolve
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
%matplotlib inline
plt.rcParams.update({"figure.dpi": 110, "font.size": 10})
# ---- constantes physiques (air, 20 °C, SI) ----
C = 343.0 # vitesse du son (m/s)
RHO = 1.204 # masse volumique (kg/m^3)
MU = 1.81e-5 # viscosité dynamique (Pa.s)
PR = 0.71 # nombre de Prandtl
GAM = 1.402 # rapport des chaleurs massiques
# géométrie de référence : col Ø2 cm, cavité Ø8 cm x 8 cm
R_COL, R_CAV, H_CAV = 0.01, 0.04, 0.08
A_H_TH = (C/(2*np.pi))*(R_COL/R_CAV)/np.sqrt(H_CAV) # constante théorique ~48.25
print(f"Constante théorique A_H^th = {A_H_TH:.2f} m^1/2/s")
Constante théorique A_H^th = 48.25 m^1/2/s
1. Théorie classique¶
Un col (longueur $L$, section $S_{col}=\pi R_{col}^2$) sur une cavité rigide $V_{cav}$ : la masse d'air du col oscille contre le ressort d'air de la cavité. La fréquence propre 1D sans perte, avec correction de bout $\Delta L$ :
$$ f_0(L) = \frac{c}{2\pi}\sqrt{\frac{S_{col}}{V_{cav}\,(L+\Delta L)}} = \frac{A_H}{\sqrt{L+\Delta L}}, \qquad A_H^{\mathrm{th}}\approx 48{,}25. $$
On cherche à estimer $\Delta L$ à partir d'un solveur 2-D, et à en discuter le statut physique.
2. Le solveur FDM axisymétrique¶
En régime harmonique, la pression complexe $p(r,z)$ obéit à l'équation de Helmholtz en coordonnées cylindriques axisymétriques ; la singularité sur l'axe de symétrie $r=0$ est levée par la règle de L'Hôpital ($\tfrac1r\partial_r p \to \partial_{rr}p$) :
$$ \partial_{rr}p + \tfrac1r\,\partial_r p + \partial_{zz}p + k^2 p = 0, \qquad k=\omega/c. $$
Conditions aux limites : Dirichlet $p=1$ à l'embouchure ($z=0$, $r\le R_{col}$), Neumann homogène sur l'axe et les parois rigides (nœuds solides réfléchis). L'assemblage de la matrice creuse est vectorisé (numpy) pour la rapidité.
In [2]:
def make_grid(w_neck=0.02, L_neck=0.04, W_cav=0.08, H_cav=0.08, h=0.002):
"""Grille structurée (r,z) et masque : 0 solide, 1 fluide, 2 embouchure (Dirichlet)."""
R_neck, R_cav = w_neck/2, W_cav/2
Nr = int(round(R_cav/h))+1; Nz = int(round((L_neck+H_cav)/h))+1
r = np.linspace(0, R_cav, Nr); z = np.linspace(0, L_neck+H_cav, Nz)
R, Z = np.meshgrid(r, z, indexing="ij")
mask = np.zeros((Nr, Nz), int)
mask[((Z < L_neck) & (R <= R_neck+1e-12)) | ((Z >= L_neck) & (R <= R_cav+1e-12))] = 1
mask[np.where(r <= R_neck+1e-12)[0], 0] = 2
return R, Z, mask, L_neck
def assemble(mask, h, k2):
"""Matrice creuse de l'opérateur de Helmholtz axisymétrique (assemblage vectorisé)."""
Nr, Nz = mask.shape; N = Nr*Nz; h2 = h*h
I = np.arange(Nr)[:,None] + Nr*np.arange(Nz)[None,:] # idx(i,j)=i+j*Nr
ii = np.broadcast_to(np.arange(Nr)[:,None], (Nr, Nz))
rows=[]; cols=[]; data=[]; b = np.zeros(N, complex)
nd = (mask != 1) # lignes Dirichlet/solide = identité
rows.append(I[nd]); cols.append(I[nd]); data.append(np.ones(int(nd.sum())))
b[I[mask == 2]] = 1.0
F = (mask == 1) # nœuds fluides intérieurs
diag = (np.where(ii == 0, -6.0/h2, -4.0/h2) + k2).astype(complex)
iis = np.where(ii == 0, 1, ii)
cl = (1 - 0.5/iis)/h2 # coeff gauche (r>0)
cr = np.where(ii > 0, (1 + 0.5/iis)/h2, 4.0/h2) # droite (axe : 4/h^2 vers i=1)
one = np.full((Nr, Nz), 1.0/h2)
diag_add = np.zeros((Nr, Nz), complex)
def shift(axis, s):
nm=np.zeros((Nr,Nz),int); nI=np.zeros((Nr,Nz),int); ok=np.zeros((Nr,Nz),bool)
if axis==0 and s>0: nm[:-1,:]=mask[1:,:]; nI[:-1,:]=I[1:,:]; ok[:-1,:]=True
if axis==0 and s<0: nm[1:,:]=mask[:-1,:]; nI[1:,:]=I[:-1,:]; ok[1:,:]=True
if axis==1 and s>0: nm[:,:-1]=mask[:,1:]; nI[:,:-1]=I[:,1:]; ok[:,:-1]=True
if axis==1 and s<0: nm[:,1:]=mask[:,:-1]; nI[:,1:]=I[:,:-1]; ok[:,1:]=True
return nm, nI, ok
def direction(coeff, src, active):
nm, nI, ok = src
real = active & ok & (nm != 0) # voisin fluide/embouchure -> hors-diag
wall = active & ~(ok & (nm != 0)) # hors-grille/solide -> repli Neumann
rows.append(I[real]); cols.append(nI[real]); data.append(coeff[real].astype(complex))
diag_add[wall] += coeff[wall]
direction(cl, shift(0,-1), F & (ii>0))
direction(cr, shift(0,+1), F)
direction(one, shift(1,-1), F)
direction(one, shift(1,+1), F)
diag = diag + diag_add
rows.append(I[F]); cols.append(I[F]); data.append(diag[F])
A = sp.coo_matrix((np.concatenate(data), (np.concatenate(rows), np.concatenate(cols))),
shape=(N, N), dtype=complex).tocsr()
return A, b, I
def solve_field(mask, h, freq):
A, b, I = assemble(mask, h, (2*np.pi*freq/C)**2)
return spsolve(A, b)[I]
def cavity_amp(mask, R, Z, L_neck, h, freq):
p = solve_field(mask, h, freq)
cav = (Z >= L_neck) & (mask == 1)
return float(np.sum(np.abs(p[cav]) * R[cav]) / np.sum(R[cav]))
print("Solveur défini.")
Solveur défini.
In [3]:
# Champ de pression à la résonance (L = 4 cm), pour visualiser le mode
R, Z, mask, Ln = make_grid(L_neck=0.04, h=0.001)
p = np.abs(solve_field(mask, 0.001, 234.6)); field = np.where(mask>0, p, np.nan)
fig, ax = plt.subplots(figsize=(3.3, 5))
pc = ax.pcolormesh(R*100, Z*100, field, shading="auto", cmap="magma")
ax.set_aspect("equal"); ax.set_xlabel("r (cm)"); ax.set_ylabel("z (cm)")
ax.set_title("Mode résonant |p| (L = 4 cm)"); fig.colorbar(pc, ax=ax, label="|p| (Pa)")
plt.show(); print(f"Grille {mask.shape} = {mask.size} nœuds")
Grille (41, 121) = 4961 nœuds
In [4]:
def f0_of(L, h, c=C):
"""Fréquence propre : repérage sur grille grossière (large) puis balayage fin + parabole."""
fpure = A_H_TH/np.sqrt(L) # borne haute (sans correction)
Rg, Zg, mg, Ln = make_grid(L_neck=L, h=0.002) # grille grossière, rapide, pour localiser
g = np.arange(0.80*fpure, 1.00*fpure, 2.0)
a = [cavity_amp(mg, Rg, Zg, Ln, 0.002, f) for f in g]
fc = g[int(np.argmax(a))] # pic approché (grille grossière)
R, Z, mask, Ln = make_grid(L_neck=L, h=h)
win = np.arange(fc-16, fc+6, 1.0) # les grilles fines résonnent un peu plus bas
aw = [cavity_amp(mask, R, Z, Ln, h, f) for f in win]
fp = win[int(np.argmax(aw))]
ff = np.arange(fp-2.5, fp+2.6, 0.25)
pf = np.array([cavity_amp(mask, R, Z, Ln, h, f) for f in ff])
i = min(max(int(np.argmax(pf)), 1), len(ff)-2); x, y = ff[i-1:i+2], pf[i-1:i+2]
A_, B_, _ = np.linalg.solve(np.column_stack((x**2, x, np.ones(3))), y)
return -B_/(2*A_)
# vérification express
print("Contrôle f0(L=4cm) : h=2mm ->", round(f0_of(0.04,0.002),2), "| h=0.5mm ->", round(f0_of(0.04,0.0005),2))
Contrôle f0(L=4cm) : h=2mm -> 234.57 | h=0.5mm -> 224.78
3. Vérification du code — solution manufacturée (MMS)¶
Sur un domaine rectangulaire lisse (sans escalier), on impose une solution analytique $p_e=\cos(a_r r)\cos(a_z z)$ (régulière sur l'axe) via un terme source $s=\mathcal{L}[p_e]$, et on mesure l'ordre de convergence. On attend l'ordre 2 (le stencil est d'ordre 2).
In [5]:
def mms_error(h):
Rd, Hd = 0.04, 0.12; ar, az = np.pi/Rd, np.pi/Hd; K = 2*np.pi*235/C
Nr=int(round(Rd/h))+1; Nz=int(round(Hd/h))+1
r=np.linspace(0,Rd,Nr); z=np.linspace(0,Hd,Nz); R,Z=np.meshgrid(r,z,indexing="ij")
pe = np.cos(ar*R)*np.cos(az*Z)
with np.errstate(divide="ignore", invalid="ignore"):
rad = np.where(R<1e-12, -ar**2*np.cos(az*Z),
-(ar/np.clip(R,1e-12,None))*np.sin(ar*R)*np.cos(az*Z))
src = (-ar**2 - az**2 + K**2)*pe + rad
N=Nr*Nz; I=np.arange(N).reshape(Nz,Nr).T; h2=h*h
jj = np.broadcast_to(np.arange(Nz)[None,:],(Nr,Nz))
bnd = (np.broadcast_to(np.arange(Nr)[:,None],(Nr,Nz))==Nr-1)|(jj==0)|(jj==Nz-1)
rows=[];cols=[];dat=[];b=np.zeros(N)
for i in range(Nr):
for j in range(Nz):
ci=I[i,j]
if bnd[i,j]: rows+=[ci];cols+=[ci];dat+=[1.0]; b[ci]=pe[i,j]; continue
if i==0:
rows+=[ci,ci,ci,ci];cols+=[I[1,j],ci,I[0,j-1],I[0,j+1]];dat+=[4/h2,-6/h2+K**2,1/h2,1/h2]
else:
rows+=[ci,ci,ci,ci,ci];cols+=[I[i-1,j],I[i+1,j],ci,I[i,j-1],I[i,j+1]]
dat+=[(1-0.5/i)/h2,(1+0.5/i)/h2,-4/h2+K**2,1/h2,1/h2]
b[ci]=src[i,j]
A=sp.coo_matrix((dat,(rows,cols)),shape=(N,N)).tocsr(); p=spsolve(A,b).reshape(Nz,Nr).T
interior=~bnd
return np.sqrt(np.mean((p[interior]-pe[interior])**2))
hs=[0.004,0.002,0.001,0.0005]; errs=[mms_error(h) for h in hs]
orders=[np.nan]+[np.log(errs[i-1]/errs[i])/np.log(2) for i in range(1,4)]
p_glob=np.polyfit(np.log(hs),np.log(errs),1)[0]
print(f"{'h (mm)':>8}{'erreur L2':>14}{'ordre':>8}")
for h,e,o in zip(hs,errs,orders):
print(f"{h*1e3:>8.2f}{e:>14.2e}{'' if o!=o else f'{o:>8.2f}'}")
print(f"\nOrdre global (régression) = {p_glob:.2f} -> CODE vérifié (ordre 2)")
h (mm) erreur L2 ordre
4.00 7.49e-03
2.00 1.81e-03 2.05
1.00 4.45e-04 2.03
0.50 1.10e-04 2.01
Ordre global (régression) = 2.03 -> CODE vérifié (ordre 2)
4. Vérification de la solution — convergence GCI (Roache)¶
Sur la vraie géométrie, on résout $f_0$ sur trois grilles ($h=2,1,0{,}5$ mm) pour cinq longueurs. Ordre observé $p$, extrapolation de Richardson $f_{h\to0}$ et incertitude GCI de la grille fine ($F_s=1{,}25$). (Calcul le plus long du notebook : ~1–2 min.)
In [6]:
Ls = [0.02, 0.03, 0.04, 0.05, 0.06]; hs3 = [0.002, 0.001, 0.0005]; Fs = 1.25
rowsG = []
for L in Ls:
f3, f2, f1 = [f0_of(L, h) for h in hs3] # grossier, moyen, fin
p = np.log(abs((f3-f2)/(f2-f1)))/np.log(2)
fext = f1 + (f1-f2)/(2**p - 1)
gci = Fs*abs((f2-f1)/f1)/(2**p - 1)
rowsG.append([L*100, f1, p, fext, gci*100])
import pandas as pd
dfG = pd.DataFrame(rowsG, columns=["L (cm)","f_fine (Hz)","p","f_extrap (Hz)","GCI (%)"]).round(
{"L (cm)":0,"f_fine (Hz)":1,"p":3,"f_extrap (Hz)":1,"GCI (%)":2})
display(dfG)
print(f"Ordre observé moyen p ≈ {dfG['p'].mean():.2f} (dégradé vers 1 par la géométrie en escalier)")
print(f"Incertitude numérique GCI ≈ {dfG['GCI (%)'].mean():.1f} %")
Lscale = np.array(Ls); f_ext = dfG["f_extrap (Hz)"].values; gci_rel = dfG["GCI (%)"].values/100
| L (cm) | f_fine (Hz) | p | f_extrap (Hz) | GCI (%) | |
|---|---|---|---|---|---|
| 0 | 2.0 | 297.0 | 0.913 | 291.9 | 2.14 |
| 1 | 3.0 | 253.6 | 0.949 | 249.6 | 1.99 |
| 2 | 4.0 | 224.8 | 0.867 | 220.6 | 2.34 |
| 3 | 5.0 | 203.9 | 0.943 | 200.6 | 2.04 |
| 4 | 6.0 | 187.9 | 0.933 | 184.7 | 2.07 |
Ordre observé moyen p ≈ 0.92 (dégradé vers 1 par la géométrie en escalier) Incertitude numérique GCI ≈ 2.1 %
In [7]:
# Contrôle du régime asymptotique : 4e grille h=0,25 mm (L=2 et 4 cm)
for L in [0.02, 0.04]:
f = [f0_of(L,h) for h in [0.002,0.001,0.0005,0.00025]]
def tri(a,b,c):
p=np.log(abs((a-b)/(b-c)))/np.log(2); return p, c+(c-b)/(2**p-1)
pA,eA = tri(f[0],f[1],f[2]); pB,eB = tri(f[1],f[2],f[3])
print(f"L={L*100:.0f}cm : p {pA:.3f}->{pB:.3f} | f_extrap {eA:.2f}->{eB:.2f} Hz "
f"(écart {abs(eA-eB):.2f} Hz -> extrapolation stable)")
L=2cm : p 0.913->0.996 | f_extrap 291.86->292.43 Hz (écart 0.57 Hz -> extrapolation stable)
L=4cm : p 0.867->1.013 | f_extrap 220.57->221.37 Hz (écart 0.80 Hz -> extrapolation stable)
5. Validation — mesures publiées (Selamet et al., JASA 1997)¶
Confrontation à des mesures publiées sur résonateurs cylindriques concentriques (col $d_c=4{,}044$ cm, $l_c=8{,}5$ cm, $V=4500$ cm³), reproduites sans paramètre ajusté. À la résonance, l'impédance d'entrée du col s'annule → la condition de Dirichlet est adéquate. (Première validation : deux configurations, $f_0$ seule.)
In [8]:
def cavity_dims(ld, V=4.5e-3):
dv = (4*V/(np.pi*ld))**(1/3); return dv, ld*dv
def f0_selamet(ld, h):
dv, lv = cavity_dims(ld)
R,Z,mask,Ln = make_grid(w_neck=0.04044, L_neck=0.085, W_cav=dv, H_cav=lv, h=h)
fpure = (C/(2*np.pi))*np.sqrt((np.pi*0.02022**2)/(4.5e-3*0.085))
g=np.arange(fpure*0.75, fpure*1.15, 2.0); a=[cavity_amp(mask,R,Z,Ln,h,f) for f in g]
fp=g[int(np.argmax(a))]; ff=np.arange(fp-3,fp+3,0.3)
pf=np.array([cavity_amp(mask,R,Z,Ln,h,f) for f in ff]); i=min(max(int(np.argmax(pf)),1),len(ff)-2)
x,y=ff[i-1:i+2],pf[i-1:i+2]; A_,B_,_=np.linalg.solve(np.column_stack((x**2,x,np.ones(3))),y)
return -B_/(2*A_)
sel=[("l/d = 1,0",1.0,91.0),("l/d = 10",10.0,72.0)]
rows=[]
for name,ld,fref in sel:
fdm=f0_selamet(ld,0.001); rows.append([name,fref,round(fdm,1),round(abs(fdm-fref)/fref*100,1)])
display(pd.DataFrame(rows, columns=["Config.","f_r publié (Hz)","FDM h=1mm (Hz)","écart (%)"]))
print("Écarts du même ordre que l'incertitude GCI (~2 %) et l'arrondi des valeurs publiées.")
| Config. | f_r publié (Hz) | FDM h=1mm (Hz) | écart (%) | |
|---|---|---|---|---|
| 0 | l/d = 1,0 | 91.0 | 93.0 | 2.2 |
| 1 | l/d = 10 | 72.0 | 72.7 | 0.9 |
Écarts du même ordre que l'incertitude GCI (~2 %) et l'arrondi des valeurs publiées.
6. Loi d'échelle (régression pondérée par l'incertitude)¶
On confronte à armes égales une loi de puissance $A\,L^b$ et le modèle physique à longueur effective $A_H/\sqrt{L+\Delta L_{\mathrm{eff}}}$, ajustés sur les fréquences extrapolées, pondérés par l'incertitude GCI.
Point méthodologique : le GCI est une borne prudente, pas un écart-type gaussien (erreurs vraisemblablement corrélées) → le $\chi^2$/AICc se lit en analyse de sensibilité. La validation leave-one-out (LOO), indépendante de ce modèle d'erreur, complète la comparaison.
In [9]:
def mpow(L,a,b): return a*L**b
def meff(L,a,dL): return a/np.sqrt(L+dL)
sig = gci_rel * f_ext # incertitude absolue (Hz)
def wfit(model, p0):
po,pc = curve_fit(model, Lscale, f_ext, p0=p0, sigma=sig, absolute_sigma=True, maxfev=30000)
chi2 = float(np.sum(((f_ext-model(Lscale,*po))/sig)**2)); dof=len(Lscale)-2
aicc = chi2 + 2*2 + 2*2*3/(len(Lscale)-2-1)
loo=[]
for i in range(len(Lscale)):
m=np.ones(len(Lscale),bool); m[i]=False
pp,_=curve_fit(model,Lscale[m],f_ext[m],p0=po,sigma=sig[m],absolute_sigma=True,maxfev=30000)
loo.append(f_ext[i]-model(Lscale[i],*pp))
return po, chi2/dof, aicc, float(np.sqrt(np.mean(np.array(loo)**2)))
poA,chiA,aicA,looA = wfit(mpow,[A_H_TH,-0.5])
poB,chiB,aicB,looB = wfit(meff,[A_H_TH,0.006])
display(pd.DataFrame({
"Modèle":["(A) puissance A·L^b","(B) longueur effective A_H/√(L+ΔL)"],
"paramètres":[f"A={poA[0]:.1f}, b={poA[1]:.3f}", f"A_H={poB[0]:.1f}, ΔL={poB[1]/R_COL:.2f} R_col"],
"χ²/dof":[f"{chiA:.2f}", f"{chiB:.3f}"], "AICc":[f"{aicA:.1f}", f"{aicB:.1f}"],
"LOO-RMSE (Hz)":[f"{looA:.2f}", f"{looB:.2f}"]}))
print(f"ΔAICc (A−B) = {aicA-aicB:.2f} -> NON décisif : les deux modèles sont compatibles dans")
print(f"l'incertitude. Le modèle physique (B) est préféré : A_H≈{A_H_TH:.1f}, ΔL={poB[1]/R_COL:.2f} R_col")
print(f"(plage de Rayleigh 0,6–0,85), et prédiction LOO {looB:.1f} vs {looA:.1f} Hz.")
Ld=np.linspace(Lscale.min()*0.9,Lscale.max()*1.1,200)
fig,(a1,a2)=plt.subplots(1,2,figsize=(9,3.5),gridspec_kw={"width_ratios":[1.5,1]})
a1.errorbar(Lscale*100,f_ext,yerr=sig,fmt="ko",capsize=3,label="FDM extrapolé ± GCI")
a1.plot(Ld*100,mpow(Ld,*poA),"r-",label=f"(A) puissance (LOO {looA:.1f})")
a1.plot(Ld*100,meff(Ld,*poB),"b--",label=f"(B) long. effective (LOO {looB:.1f})")
a1.set_xlabel("L (cm)");a1.set_ylabel("f₀ (Hz)");a1.legend(fontsize=8);a1.grid(ls=":")
for mdl,po,c,lab in [(mpow,poA,"r","(A)"),(meff,poB,"b","(B)")]:
a2.errorbar(Lscale*100,f_ext-mdl(Lscale,*po),yerr=sig,fmt="o",color=c,capsize=3,label=lab)
a2.axhline(0,color="k",lw=.8);a2.set_xlabel("L (cm)");a2.set_ylabel("résidu (Hz)")
a2.legend(fontsize=8);a2.grid(ls=":");plt.tight_layout();plt.show()
| Modèle | paramètres | χ²/dof | AICc | LOO-RMSE (Hz) | |
|---|---|---|---|---|---|
| 0 | (A) puissance A·L^b | A=57.4, b=-0.417 | 0.09 | 10.3 | 2.97 |
| 1 | (B) longueur effective A_H/√(L+ΔL) | A_H=47.7, ΔL=0.67 R_col | 0.005 | 10.0 | 0.58 |
ΔAICc (A−B) = 0.26 -> NON décisif : les deux modèles sont compatibles dans l'incertitude. Le modèle physique (B) est préféré : A_H≈48.3, ΔL=0.67 R_col (plage de Rayleigh 0,6–0,85), et prédiction LOO 0.6 vs 3.0 Hz.
7. Décomposition de la correction de bout (impédance de rayonnement)¶
L'embouchure Dirichlet ne contient pas la charge de rayonnement extérieure : $\Delta L_{\mathrm{eff}}$ est une correction effective. On construit une seconde configuration où l'embouchure porte l'impédance de rayonnement d'un piston bafflé, $Z_{\mathrm{rad}}=\rho c[\tfrac12(kR_{col})^2+i\tfrac{8}{3\pi}kR_{col}]$, dont la réactance encode $\Delta L_{\mathrm{ext}}=8R_{col}/(3\pi)\approx0{,}849\,R_{col}$. La différence des longueurs effectives isole intérieur/extérieur.
⚠️ La réactance imposée contient $8R_{col}/(3\pi)$ : c'est un test de cohérence de la chaîne (condition de Robin bien implémentée, régression fidèle), pas un calcul ab initio.
In [10]:
def solve_config(L, freq, radiation, h=0.001):
"""Réponse de cavité avec forçage volumique en cavité ; l'embouchure est soit Dirichlet
p=0 (intérieur seul), soit une impédance de rayonnement de piston bafflé (Robin)."""
R, Z, mask, Ln = make_grid(L_neck=L, h=h); Nr, Nz = mask.shape
omega = 2*np.pi*freq; k2 = (omega/C)**2; inv_h2 = 1/h**2; N = Nr*Nz
idx = lambda i, j: i + j*Nr
rows=[]; cols=[]; data=[]; b=np.zeros(N, complex)
zc = L + 0.5*H_CAV; F = 1e4*np.exp(-((R)**2 + (Z-zc)**2)/(2*0.01**2)) # source volumique
fluid = (mask>0); mouth = (mask==2)
ka = omega/C*R_COL; Zr = RHO*C*(0.5*ka**2 + 1j*(8/(3*np.pi))*ka); alpha = 1j*omega*RHO/Zr
for i in range(Nr):
for j in range(Nz):
ci = idx(i, j)
if not fluid[i, j]: rows+=[ci]; cols+=[ci]; data+=[1.0]; continue
if mouth[i, j]:
if not radiation: rows+=[ci]; cols+=[ci]; data+=[1.0]; b[ci]=0.0; continue
if i==0:
diag=-6*inv_h2+k2-2*alpha/h
rows+=[ci]; cols+=[idx(1,j)]; data+=[4*inv_h2]
rows+=[ci]; cols+=[idx(0,j+1)]; data+=[2*inv_h2]
else:
diag=-4*inv_h2+k2-2*alpha/h; lc=inv_h2*(1-0.5/i); rc=inv_h2*(1+0.5/i)
if not fluid[i-1,j]: diag+=lc
else: rows+=[ci]; cols+=[idx(i-1,j)]; data+=[lc]
if i+1>=Nr or not fluid[i+1,j]: diag+=rc
else: rows+=[ci]; cols+=[idx(i+1,j)]; data+=[rc]
rows+=[ci]; cols+=[idx(i,j+1)]; data+=[2*inv_h2]
rows+=[ci]; cols+=[ci]; data+=[diag]; b[ci]=-F[i,j]; continue
if i==0:
diag=-6*inv_h2+k2; rows+=[ci]; cols+=[idx(1,j)]; data+=[4*inv_h2]
for nj in (j-1,j+1):
if nj<0 or nj>=Nz or not fluid[0,nj]: diag+=inv_h2
else: rows+=[ci]; cols+=[idx(0,nj)]; data+=[inv_h2]
else:
diag=-4*inv_h2+k2; lc=inv_h2*(1-0.5/i); rc=inv_h2*(1+0.5/i)
if not fluid[i-1,j]: diag+=lc
else: rows+=[ci]; cols+=[idx(i-1,j)]; data+=[lc]
if i+1>=Nr or not fluid[i+1,j]: diag+=rc
else: rows+=[ci]; cols+=[idx(i+1,j)]; data+=[rc]
for nj in (j-1,j+1):
if nj<0 or nj>=Nz or not fluid[i,nj]: diag+=inv_h2
else: rows+=[ci]; cols+=[idx(i,nj)]; data+=[inv_h2]
rows+=[ci]; cols+=[ci]; data+=[diag]; b[ci]=-F[i,j]
A = sp.coo_matrix((data,(rows,cols)), shape=(N,N), dtype=complex).tocsr()
p = spsolve(A, b).reshape((Nz,Nr)).T
cav = (Z>=L) & (mask==1); return float(np.sum(np.abs(p[cav])*R[cav])/np.sum(R[cav]))
def f0_config(L, radiation):
fc = A_H_TH/np.sqrt(L + (0.007 if not radiation else 0.015))
g = np.arange(fc-18, fc+8, 2.0); a=[solve_config(L,f,radiation) for f in g]; fp=g[int(np.argmax(a))]
ff = np.arange(fp-2.5, fp+2.6, 0.25); pf=np.array([solve_config(L,f,radiation) for f in ff])
i=min(max(int(np.argmax(pf)),1),len(ff)-2); x,y=ff[i-1:i+2],pf[i-1:i+2]
A_,B_,_=np.linalg.solve(np.column_stack((x**2,x,np.ones(3))),y); return -B_/(2*A_)
Lr = np.array([0.02,0.03,0.04,0.05,0.06])
f_dir = np.array([f0_config(L,False) for L in Lr]); f_rad = np.array([f0_config(L,True) for L in Lr])
pd_,_ = curve_fit(meff, Lr, f_dir, p0=[48,0.007]); pr_,_ = curve_fit(meff, Lr, f_rad, p0=[48,0.015])
shift = pr_[1]-pd_[1]; ana = 8*R_COL/(3*np.pi)
display(pd.DataFrame({"Embouchure":["Dirichlet (intérieur seul)","Rayonnement (total)","Décalage (extérieur)"],
"ΔL / R_col":[f"{pd_[1]/R_COL:.3f}", f"{pr_[1]/R_COL:.3f}", f"{shift/R_COL:.3f}"]}))
print(f"Décalage {shift/R_COL:.3f} R_col vs analytique 8/(3π) = {ana/R_COL:.3f} R_col -> écart {abs(shift-ana)/ana*100:.1f} %")
| Embouchure | ΔL / R_col | |
|---|---|---|
| 0 | Dirichlet (intérieur seul) | 0.676 |
| 1 | Rayonnement (total) | 1.499 |
| 2 | Décalage (extérieur) | 0.823 |
Décalage 0.823 R_col vs analytique 8/(3π) = 0.849 R_col -> écart 3.1 %
8. Pertes viscothermiques et facteur de qualité $Q$¶
Le solveur fréquentiel est sans perte. Deux mécanismes d'amortissement : l'absorption volumique de Stokes (négligeable, $Q\sim10^6$) et la friction viscothermique de couche limite dans le col (dominante). Bilan d'énergie (résistance de surface $R_s=\sqrt{\rho\mu\omega/2}\,[1+(\gamma-1)/\sqrt{Pr}]$) :
$$ Q_{\mathrm{visc}}=\frac{R_{col}}{F_t\,\delta_v}, \quad \delta_v=\sqrt{\tfrac{2\mu}{\rho\omega_0}}, \quad F_t=1+\tfrac{\gamma-1}{\sqrt{Pr}}\approx1{,}48. $$
On l'illustre par un ring-down amorti (schéma leapfrog explicite) : un amortissement localisé dans le col, dont l'amplitude est calibrée pour reproduire $Q_{\mathrm{visc}}$ (contrôle de cohérence).
In [11]:
f0 = 234.6; w0 = 2*np.pi*f0
Ft = 1 + (GAM-1)/np.sqrt(PR)
dv = np.sqrt(2*MU/(RHO*w0)); Q_visc = R_COL/(Ft*dv)
delta_diff = (1/RHO)*((4/3)*MU + 0.6*MU + (GAM-1)*0.0257/1005)
Q_bulk = w0/(delta_diff*w0**2/C**2)
print(f"δ_v = {dv*1e6:.0f} µm | Q_visc (couche limite col) ≈ {Q_visc:.1f} | Q_bulk (Stokes) ≈ {Q_bulk:.1e} (négligeable)")
# --- ring-down transitoire (leapfrog) avec amortissement localisé dans le col ---
R,Z,mask,Ln = make_grid(L_neck=0.04, h=0.002); Nr,Nz=mask.shape
dt = 0.85*0.002/(C*np.sqrt(2)); Nt=int(0.05/dt); h2=0.002**2; c2dt2=(C*dt)**2
r=R; inv_r=np.where(r<1e-7,0,1/np.clip(r,1e-7,None)); neck=(Z<0.04)&(mask>0)
mL=np.zeros_like(mask,bool);mL[1:]=(mask[:-1]==0); mR=np.zeros_like(mask,bool);mR[:-1]=(mask[1:]==0)
mB=np.zeros_like(mask,bool);mB[:,1:]=(mask[:,:-1]==0); mT=np.zeros_like(mask,bool);mT[:,:-1]=(mask[:,1:]==0)
def ring(gamma_amp):
g=np.where(neck,gamma_amp,0.0); pp=np.zeros((Nr,Nz)); pc=np.zeros((Nr,Nz)); sig=[]
den=1+g*dt/2
for n in range(Nt):
t=n*dt; pad=np.pad(pc,1,mode="edge")
L_=np.where(mL,pc,pad[0:Nr,1:Nz+1]); Rr=np.where(mR,pc,pad[2:Nr+2,1:Nz+1])
B_=np.where(mB,pc,pad[1:Nr+1,0:Nz]); T_=np.where(mT,pc,pad[1:Nr+1,2:Nz+2])
d2r=(Rr-2*pc+L_)/h2; dr=(Rr-L_)/(2*0.002)
lap=d2r+np.where(r<1e-7,d2r,dr*inv_r)+(T_-2*pc+B_)/h2
pn=(2*pc-pp+c2dt2*lap+(g*dt/2)*pp)/den
pn[mask==2]=np.sin(w0*t)*np.exp(-((t-0.003)**2)/(2*0.001**2)); pn[mask==0]=0
cav=(Z>=0.04)&(mask==1); sig.append(float(np.sum(pn[cav]*R[cav])/np.sum(R[cav])))
pp,pc=pc,pn
return np.array(sig)
def measure_Q(sig):
from numpy import polyfit
t=np.arange(len(sig))*dt; m=t>=0.008
from scipy.signal import hilbert
env=np.abs(hilbert(sig[m])); tt=t[m]; good=env>0.05*env.max()
b=-polyfit(tt[good],np.log(env[good]),1)[0]; return w0/(2*b)
# calibrage par point fixe (Q ~ K/amp) pour atteindre Q_visc
amp = 300.0
for _ in range(5):
Qm = measure_Q(ring(amp))
if abs(Qm - Q_visc)/Q_visc < 0.03: break
amp = Qm*amp/Q_visc
sig_visc = ring(amp); sig_none = ring(0.0); Qm = measure_Q(sig_visc)
t_ms=np.arange(len(sig_visc))*dt*1e3
fig,ax=plt.subplots(figsize=(7.5,3.2))
ax.plot(t_ms,sig_none,color="0.6",lw=1,label="sans perte")
ax.plot(t_ms,sig_visc,"b-",lw=1.2,label=f"viscothermique col (Q={Qm:.0f})")
ax.set_xlabel("t (ms)");ax.set_ylabel("pression moyenne cavité (Pa)")
ax.set_title("Ring-down : décroissance amortie");ax.legend(fontsize=8);ax.grid(ls=":")
plt.tight_layout();plt.show()
print(f"Q analytique {Q_visc:.1f} vs Q mesuré (décrément) {Qm:.1f} — contrôle de cohérence "
f"(écart {abs(Qm-Q_visc)/Q_visc*100:.0f} %)")
δ_v = 143 µm | Q_visc (couche limite col) ≈ 47.4 | Q_bulk (Stokes) ≈ 2.1e+06 (négligeable)
Q analytique 47.4 vs Q mesuré (décrément) 47.4 — contrôle de cohérence (écart 0 %)
9. Conclusion¶
- V&V hiérarchisée : code (MMS, ordre 2,03) → solution (GCI ≈2 %, 4e grille de contrôle) → modèle (mesures publiées, écarts 0,6 % et 2,2 %).
- Loi d'échelle : loi de puissance et longueur effective indiscernables dans l'incertitude ($\Delta$AICc non décisif) ; le modèle physique est préféré (constante théorique, $\Delta L_{\mathrm{eff}}\approx0{,}66\,R_{col}$, LOO 10× meilleure). L'exposant $b\approx-0{,}42$ est un artefact de plage restreinte.
- Correction de bout décomposée intérieur/extérieur (accord à ~3 % avec $8R_{col}/(3\pi)$, test de cohérence).
- Pertes dominées par la couche limite du col ($Q\approx47$) ; absorption volumique négligeable.
Perspectives : domaine extérieur maillé (Sommerfeld/PML) pour un rayonnement ab initio ; validation élargie (plusieurs rapports col/cavité, réponse complète) ; campagne expérimentale propre.
Notebook autonome — aucune dépendance locale, uniquement numpy / scipy / matplotlib.