% Comparison of models by Delany-Bazley 
% and Champoux-Allard to compute the sound
% absorption coefficient of a material sample
% backed by an impervious rigid wall.
%
% M. E. Delany and E. N. Bazley, 
% Acoustical properties of fibrous absorbent materials,
% Applied Acoustics (3), 1970, pp. 105-116
%
% J.-F. Allard and Y. Champoux, 
% New empirical equations for sound propagation in rigid frame
%   fibrous materials,
% J. Acoust. Soc. Am. 91(6), 1992, pp. 3346-3353
%
% Copyleft 2006 luc.jaouen@matelys.com
% cf. APMR on the web,
% PropagationModels/MotionlessSkeleton/DelanyBazleyModel.html
% for more information

close all
clear all

f = [200:5:8000];   % [s-1] frequency range
omega = 2*pi*f;     % [s-1] angular frequency

rho_0 = 1.213;      % [Kg.m-3] density at rest of air at 18C, 1atm
c_0   = 342.2;      % [m.s-1] speed of sound in air at 18C, 1atm
P_0   = 1.0132e+05; % [N.m-2] atmospheric pressure at 18C, 1atm
sigma = 10000       % [N.s.m-4] static air flow resistivity of material
h     = 0.05        % [m] thickness of material sample


%%%%%
%%%%% Compute variable X and print frequency
%%%%% limits of validity for the two models 
%%%%%
X = f/sigma;
f_min = 0.01*sigma
f_max = 1.00*sigma


%%%%%
%%%%% Delany and Bazley model
%%%%% (NB: gamma = alpha + j beta = j k )
%%%%%
Z_DB70 = rho_0*c_0*( 1 + 9.08*(X*1000).^(-0.75) ...
                     - i*11.9*(X*1000).^(-0.73) ); 

k_DB70 = omega/c_0 .* (-i) .* ( 10.3*(X*1000).^(-0.59) ...
                                + i* ( 1 + 10.8*(X*1000).^(-0.70) ) );

K_DB70 = Z_DB70.*omega./k_DB70;
rho_DB70 = k_DB70.*Z_DB70./omega;


%%%%%
%%%%% Expressions of Delany and Bazley model given by
%%%%% Allard and Champoux
%%%%% 

Z_AC92 = rho_0*c_0*( 1 + 0.0571*(rho_0*X).^(-0.754) ...
                     - i*0.087*(rho_0*X).^(-0.732) ); 

k_AC92 = omega/c_0 .* ( 1+0.0978*(rho_0*X).^(-0.700) ...
                        - i*0.189*(rho_0*X).^(-0.595) );

K_AC92 = Z_AC92.*omega./k_AC92;
rho_AC92 = k_AC92.*Z_AC92./omega;


%%%%%
%%%%% Compute sound absorption using the two models 
%%%%% for a sample of thickness d backed by a rigid 
%%%%% and impervious wall under at room temperature
%%%%% and pressure conditions 
%%%%%

Z = -j.*Z_DB70./tan(k_DB70*h);
alpha_DB70 = 1 - ( abs( (Z-rho_0*c_0)./(Z+rho_0*c_0) ) ).^2;


Z = -j.*Z_AC92./tan(k_AC92*h);
alpha_AC92 = 1 - ( abs( (Z-rho_0*c_0)./(Z+rho_0*c_0) ) ).^2;


%%%%%
%%%%% Compare results
%%%%%

figure(1)
set(gca,'FontSize',16)
plot(f,real(rho_DB70)/rho_0,'k-','LineWidth',2)
hold on
plot(f,imag(rho_DB70)/rho_0,'b-','LineWidth',2)
plot(f,real(rho_AC92)/rho_0,'r--','LineWidth',2)
plot(f,imag(rho_AC92)/rho_0,'m--','LineWidth',2)
xlabel('Frequency (Hz)')
ylabel('Normalized dynamic density')
legend('DB70 (Re)','DB70 (Im)','AC92 (Re)','AC92 (Im)',4)

figure(2)
set(gca,'FontSize',16)
plot(f,real(K_DB70)/P_0,'k-','LineWidth',2);
hold on
plot(f,imag(K_DB70)/P_0,'b-','LineWidth',2);
plot(f,real(K_AC92)/P_0,'r--','LineWidth',2);
plot(f,imag(K_AC92)/P_0,'m--','LineWidth',2);
xlabel('Frequency (Hz)')
ylabel('Normalized dynamic bulk modulus')
legend('DB70 (Re)','DB70 (Im)','AC92 (Re)','AC92 (Im)',4)

figure(3)
set(gca,'FontSize',16)
plot(f,alpha_DB70,'k-','LineWidth',2);
hold on
plot(f,alpha_AC92,'r--','LineWidth',2);
xlabel('Frequency (Hz)')
ylabel('Sound absorption coefficient')
legend('DB70','AC92',4)