In 1956, M. A. Biot wrote two articles [Bio56a, Bio56b] where he proposes a theoretical formulism to describe the elastic wave propagation inside an isotropic porous medium saturated with a fluid. One of the key result of this work is the existence of three waves : two compressional waves and one shear wave.

There are already numerous books, theses or articles introducing the Biot's theory for acoustic materials. This page focuses on an essential point and two anecdotic comments which are [almost] never discussed.
An explorable explanation and a commented Matlab/GNU Octave script are added to offer a different "perspective" to the reader.

The fastest compressional wave is usally called the $P_1$ wave. The slowest compressional wave is usually called the $P_2$ wave. The shear wave is called, as for isotropic elastic media, the $S$ wave.

The long wavelength (compare to the pore dimensions) and the small deformation hypotheses ensure that, when Biot's theory can be applied, the fluid properties (dynamic mass density & dynamic bulk modulus) can be considered as equivalent to their values when the solid phase is not in motion. In other words, the deformation of the pore skeleton does not significantly affect the fluid properties.

Biot's work, first developed for geophysics & oil prospection, led to controversy during a time, for multiple reasons.
First due to the technique used to derivate the wave equations : a form of the Lagrangian was postulated (before the use of Hamilton's principle). Indeed, the homogenization theory has been developed later on. Second, the main results of Biot's work is the prediction of a slow compressional wave $P_2$ which was difficult to observe at low (infrasounds) or high frequencies (ultrasounds) until the ultrasonic experimental evidence by Plona in 1980 [Plo80] who observed travel times of both compressional waves consistent with Biot’s predictions from sintered glass beads in water.

It is interesting to note here that acousticians were observing, since long, two compressional waves (while not identified as primary and secondary compressional waves) when measuring the normal surface impedance, for plane waves, of porous materials. As an example, the first two sentences of Leo Beranek's 1947 publication [Ber47] read : "Studies on rigid acoustical tiles and soft blankets are described in this paper. It is shown that two waves travel through the material -- one primarily airborne and the other primarily structure-borne" and at page 564 : "The presence of resonances of the disks of material due to the clamping effect of the edges is clearly seen".

Explorable absorption example

The application below let you play with the parameters of a porous material (modeled as a Biot-JCAL or a JCAL material with a rigid skeleton) and see the influence of each parameter on the diffuse sound absorption of this material backed by a rigid & impervious wall.

The JCAL curve (■) is obtained while considering the material has a rigid and motionless skeleton (i.e. we can consider the curve as the result of the only $P_2$ wave). The Biot-JCAL curve (■) is obtained considering the $P_1$, $P_2$ and $S$ waves.

Absorption of a porous material (isotropic Biot-JCAL model vs. JCAL model)

thickness : mm,
porosity : %,
resistivity : N.s.m^-4,
tortuosity : ,
viscous charac. length : microns,
thermal charac. length : microns,
static th. permeability : × 10^-10 m²

mass density : kg m^-3
Young's modulus : Pa
Poisson's ratio :
loss factor :

Sound abs. coefficient vs. frequency (Hz):
■ rigid skeleton (JCAL model)
■ elastic skeleton (Biot-JCAL model)

The Biot-JCAL simulation is obtained using the formulation introduced in [DBJ12] (for application to double-porosity materials) and further detailled and applied to multiple models (including perforated plates) in [BJ13]. The Transfer Matrix Method (TMM) for acoustical materials is used to compute these diffuse sound field simulations, as implemented in [BLA95] for the JCAL model and [DBGP13] for the Biot-JCAL one.

One implementation of the TMM should have been sufficient. The 2 implementations are only used to get a reduction of the computation operations and to offer an almost real-time experience even on mobile phones.

Matlab/GNU Octave script explained

BiotJCAL.m is a Matlab/GNU Octave script to compute the sound absorption, for plane waves at normal incidence (it thus cannot be compared to the diffuse sound field interactive computation above) of a material sample modeled using Biot's theory together with a JCAL model.
To check the script we compare our result with one (SolutionFromAlphaCell.rok) from a Transfer Matrix Method (TMM) computation. Both results should be equivalent.

Remember this script is only valid for plane waves at normal incidence. If you want to obtain the sound absorption for a diffuse sound field, try to implement a Transfer Matrix Method (TMM) rather than starting from this script.

First, a description of the script & references used.
Note that the time convention is $e^{j\omega t}$

% Compute the sound absorption coefficient of a porous material sample
% for plane waves, at normal incidence (i.e. impedance tube conditions).
% The Biot-Johnson-Champoux-Allard-Lafarge (Biot-JCAL) model 
% is used to compute the behavior of the material sample.
%
% No shear wave (S) is accounted here due to the conditions.
%
% This script is written using a e^{+j omega t} time convention.
%
% It's not forbidden to keep the references reported in this work
% and to cite this work if you use it (there's no shame to let others
% know you built nice stuff based on the works of other nice people).
%
% Copyleft 2012 luc.jaouen@matelys.com
% cf. APMR on the web:
% apmr.matelys.com/PropagationModels/BiotTheory/
% for more information.

clear all
close all

%%%%%
%%%%% Geometry
%%%%%

Self explanatory.

h_sample = 57e-3; % [m] sample thickness

%%%%%
%%%%% Parameters of the 
%%%%% Johnson-Champoux-Allard-Lafarge (JCAL) model
%%%%%

We now define the 6 parameters of the JCAL model
The elastic parameters of the solid phase (Biot's part of the model) are introduced in a second step.

sigma        = 40000;   % [N.s.m-4] static air flow resistivity of material
phi          = 0.94;    % open porosity
alpha_infty  = 1.06;    % high frequency limit of the tortuosity
lambda       = 56e-06;  % [m] viscous characteristic length
lambda_prime = 140e-06; % [m] thermal characteristic length
k_0_prime    = 23e-10;  % [m2] static thermal permeability

%%%%%
%%%%% Elastic parameters to get Biot-JCAL model
%%%%%

rho1         = 80;      % mass density of the porous material [kg.m-3]
young        = 380000;  % Young's modulus (Pa = N.m-2]
poisson      = 0.0;     % Poisson's ratio
eta_s        = 0.2;     % loss factor 

%%%%%
%%%%% Frequency 'loop'
%%%%%

We define a frequency spectrum.
f is obviously the frequency, omega the pulsation.

f = [20:10:4000];
omega = 2*pi*f;

%%%%%
%%%%% Atmospheric conditions / constants
%%%%%

T0 = 20;      % Celsius degrees
P0 = 101325;  % Pa

%%%%%
%%%%% Do not modify anything below this line
%%%%% unless you understand what you're doing
%%%%%
% Compute air parameters related to room conditions
% (acknowledgment goes to Dominic Pilon, 
% dominic.pilon@metafoam.com, for gathering almost all these
% expressions from the references below)
% 
% References:
% Lide, D. R. and Kehiaian H. V.,
% CRC. Handbook of Thermophysical and Thermochemical Data,
% CRC. Press Inc, 1994
%
% Touloukian, Y. S. and Makita, T.,
% Specific Heat - Non metallic Liquids and gases,
% The TPRC Data Series, Volume 6, IFI/PLENUM, 1970
% 
% Pierce, A. D.,
% Acoustics, An Introduction to Its Physical Principles and Applications
% Acoustical Society of America, 2nd edition, 1989

% Convert temperature from Celsius to Kelvin
T     = 273.16+T0;
  
% Universal gas constant (J.K-1.kg-1) [also 8.314 J.mol-1.K-1]
R     = 287.031;            

% Specific heat at constant pressure (J.kg-1.K-1; 260 K < T < 600 K
Cp    = 4168.8*(0.249679-7.55179e-5*T+1.69194e-7*T^2-6.46128e-11*T^3);  
  
% Specific heat at constant volume (J.kg-1.K-1; 260 K < T < 600 K
Cv    = Cp-R;               

% Dynamic viscosity (N.s.m-2; 100 K < T < 600 K
eta   = 7.72488e-8*T-5.95238e-11*T^2+2.71368e-14*T^3;  

% Ratio of specific heats
gamma = Cp/Cv;

% Density of air (kg.m-3)
rho0 = P0/(R*T);

% Velocity of sound (m.s^-1)
c0   = sqrt(gamma*R*T);

% Thermal conductivity  (W.m-1.K-1) - cf A. D. Pierce p 513
kappa = 2.624e-02 * ( (T/300)^(3/2) * (300+245.4*exp(-27.6/300))/(T+245.4*exp(-27.6/T)) );

% Prandtl's number
Pr    = eta*Cp/kappa;

In this part, we compute the dynamic mass density $\rho_{eq}$ and the dynamic bulk modulus $K_{eq}$ of the material, for each frequency, following the JCAL model. An intermediate variable $G_j$ is used to compute $\rho_{eq}$. \[ \rho_{eq}(\omega) = \frac{ \alpha_{\infty} \rho_{0} }{ \phi } \Bigg[ 1 + \frac{ \sigma \phi }{ j \omega \rho_{0} \alpha_{\infty}} \sqrt{ 1 + j\frac{ 4 \alpha_{\infty}^{2} \eta \rho_{0} \omega }{ \sigma^{2} \Lambda^{2} \phi^{2}} }\Bigg] \] \[ K_{eq}(\omega) = \displaystyle{\frac{\gamma P_{0}/\phi} {\gamma - (\gamma-1) \left[ 1-j\displaystyle{\frac{\phi\kappa} {k'_0 C_p\rho_{0}\omega}} \sqrt{ 1+j \displaystyle{\frac{4{k'}_{0}^{2} C_p\rho_{0}\omega} {\kappa\Lambda'^2\phi^2}} } \, \right]^{-1} }} \] $K_e$ is the bulk modulus of the fluid phase that we will need later on.

%%%%%
%%%%% COMPUTE THE DYNAMIC MASS DENSITY AND BULK MODULUS FOR THE MATERIAL
%%%%% JCAL MODEL
%%%%%

Gj = sqrt( 1 + 4*j*alpha_infty^2*eta*rho0*omega / ...
             (sigma^2*lambda^2*phi^2)          );

rho_eq = (1/phi) * alpha_infty * rho0 * ( 1 + sigma*phi*Gj ./ ... 
                         (j*omega*rho0*alpha_infty) );

K_eq = (1/phi) * gamma*P0 ./ ...
       ( gamma - (gamma-1) ./ ...
       ( 1 -  j*phi*kappa./(k_0_prime*Cp*omega*rho0) .* ...
       sqrt( 1+j*4*rho0*omega*Cp*k_0_prime^2/(kappa*lambda_prime^2*phi^2) ) ) );

K_e = K_eq*phi;

We then compute the mass densities and lastic moduli needed in Biot's theory.
We use [DBJ12, BJ13] formalism to express these quantities from $\rho_{eq}$ and $K_{eq}$ (or $K_e$).
In particular, for the mass densities we have : \[ \begin{align} \rho_{22} &= \phi^2\rho_{eq} \\ \rho_{11} &= \phi^2\rho_{eq} - \phi\rho_0 + \rho_1 \\ \rho_{12} &= \phi\rho_0 - \phi^2\rho_{eq} \end{align} \]

%%%%%
%%%%% COMPUTE EQUIVALENT MASS DENSITIES AND BIOT PARAMETERS
%%%%%
% Using the formalism introduced in 
% O. Dazel, F.-X. Becot, L. Jaouen, 
% Biot effects for sound absorbing double porosity materials, 
% Acta Acustica United with Acustica 98(4), pp 567-576, 2012.
% and further detailed in
% F.-X. Becot, L. Jaouen, 
% An alternative Biot's formulation for dissipative porous media with skeleton deformation, 
% J. Acoust. Soc. Am. 134 (6), pp. 4801-4807, 2013.
%
% Equation numbers correspond to those from
% - J.-F. Allard, Propagation of sound in porous media, first ed., 1993
% - or second ed., 2009, with N. Atalla

rho_22_tilde = rho_eq*phi^2;
rho_11_tilde = rho_eq*phi^2-phi*rho0+rho1;
rho_12_tilde = phi*rho0-rho_eq*phi^2;

N = (young*(1+i*eta_s))/(2*(1+poisson));  
KB = (2/3)*N*(1+poisson)/(1-2*poisson);      % eq. 6.29

rho_s_tilde = rho_11_tilde-((rho_12_tilde.^2)./rho_22_tilde);

P_tilde = (((1-phi)^2)/phi).*K_e+KB+(4/3)*N; % eq. 6.28
Q_tilde = (1-phi)*K_e;                       % eq. 6.27
R_tilde = phi*K_e;                           % eq. 6.26

%%%%%
%%%%% COMPUTE COMPLEX WAVE NUMBERS
%%%%% 
delta = (((P_tilde.*rho_22_tilde)+(R_tilde.*rho_11_tilde)-(2.*Q_tilde.*rho_12_tilde)).^2) ...  % eq. 6.71
      -(4.*(P_tilde.*R_tilde-Q_tilde.^2).*(rho_11_tilde.*rho_22_tilde-rho_12_tilde.^2));

sdelta = (sqrt(delta));

% Take square root of delta having a positive real part
if real(sdelta) < 0
    sdelta = - sdelta;
end

delta_1 = sqrt((omega.^2./(2.*(P_tilde.*R_tilde-Q_tilde.^2)))... % eq. 6.69
          .*(P_tilde.*rho_22_tilde+R_tilde.*rho_11_tilde-2.*Q_tilde.*rho_12_tilde+(sdelta)));

delta_2 = sqrt((omega.^2./(2.*(P_tilde.*R_tilde-Q_tilde.^2)))... % eq. 6.70
          .*(P_tilde.*rho_22_tilde+R_tilde.*rho_11_tilde-2.*Q_tilde.*rho_12_tilde-(sdelta)));

mu_1 = ((P_tilde.*delta_1.^2)-(rho_11_tilde.*omega.^2))./((rho_12_tilde.*omega.^2)-(Q_tilde.*delta_1.^2)); % eq. 6.73
mu_2 = ((P_tilde.*delta_2.^2)-(rho_11_tilde.*omega.^2))./((rho_12_tilde.*omega.^2)-(Q_tilde.*delta_2.^2)); % eq. 6.73

%%%%%
%%%%% CHARACTERISTIC IMPEDANCES OF WAVES FOR SOLID AND FLUID PHASES
%%%%% 

Z1f = (R_tilde+Q_tilde./mu_1).*delta_1./(phi*omega); % eq. 6.75
Z2f = (R_tilde+Q_tilde./mu_2).*delta_2./(phi*omega); % eq. 6.76
Z1s = (P_tilde+Q_tilde.*mu_1).*delta_1./omega;       % eq. 6.78
Z2s = (P_tilde+Q_tilde.*mu_2).*delta_2./omega;       % eq. 6.79

%%%%%
%%%%% COMPUTATION OF THE SURFACE IMPEDANCE
%%%%%

D = (1-phi+phi*mu_2).*(Z1s-(1-phi).*Z1f.*mu_1).*tan(delta_2*h_sample) + ...
    (1-phi+phi*mu_1).*(Z2f.*mu_2.*(1-phi)-Z2s).*tan(delta_1*h_sample); % eq. 6.109

Zs = -1j*(Z1s.*Z2f.*mu_2-Z2s.*Z1f.*mu_1)./D; % eq. 6.108

\[ \alpha = 1 - |R|^2 = 1 - \left|\frac{Z_s - Z_0}{Z_s + Z_0}\right|^2 \]

%%%%%
%%%%% COMPUTE SOUND ABSORPTION COEFFICIENT
%%%%%
alpha  = 1 - ( abs( (Zs-rho0*c0)./(Zs+rho0*c0) ) ).^2;

To check the script we compare our result with one from a Transfer Matrix Method (TMM) computation. Both results should be equivalent.

alphacell = load('SolutionFromAlphaCell.rok');

%%%%%
%%%%% PLOT THE RESULT
%%%%%

figure(1)
set(gca,'FontSize',16)
hold on
plot(f,alpha,'r-','Linewidth',2)
plot(alphacell(:,1),alphacell(:,2),'y--','Linewidth',2)
xlabel('Frequency (Hz)')
title('Normal sound incidence (Biot-JCAL model)')
set(gca,'Ylim',[0 1])