% 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
%%%%%

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

%%%%%
%%%%% The material
%%%%% (Johnson-Champoux-Allard-Lafarge (JCAL) model)
%%%%%

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
%%%%%

f = [20:5: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;


%%%%%
%%%%% 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;

%%%%%
%%%%% 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;                % eq. 6.44
rho_11_tilde = rho_eq*phi^2-phi*rho0+rho1; % eq. 6.40
rho_12_tilde = phi*rho0-rho_eq*phi^2;       % eq. 6.42


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_2 = 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_1 = 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


%delta_3 = omega.*sqrt(rho_s_tilde./N);
%mu_3 = -rho_12_tilde./rho_22_tilde;

%%%%%
%%%%% 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

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

%%%%%
%%%%% COMPUTE SOUND ABSORPTION COEFFICIENT
%%%%%

alpha  = 1 - ( abs( (Zs-rho0*c0)./(Zs+rho0*c0) ) ).^2;

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])

