Post
by **Jim Kirby** » Sat Dec 29, 2012 12:17 am

I couldn't find a file extension that the board would let me upload (what do you use for ASCII text files?

Here is the code inline.

Jim

%--------------------------------------------------------------------------

%

% gore_gillet_4dof.m

%

% Implements the 4 DOF model described in section 2.4.1 in Gore and Gillet,

% etc.

%

%

% copyright Jim Kirby, 11/27/12 - 12/28/12, free for use by anyone.

%

%--------------------------------------------------------------------------

%--------------------------------------------------------------------------

%

% Physical constants(Table 2.4-1, pg 2-37.

%

%--------------------------------------------------------------------------

mt=0.07; % mass of top

mb=0.094; % mass of back

ms=0.8; % mass of sides

ma=0.000591; % sound hole air mass

At=0.0430; % area of top

Ab=0.04; % area of back

Aa=0.0063; % soundhole area

As=0.025; % area of sides

Rt=3.0; % damping coefficient, top

Rb=7.0; % damping coefficient, back

Rs=10.0; % damping coefficient, sides

Ra=0.014; % damping coefficient, air

Kt=68000; % equivalent top plate stiffness

Kb=130000; % equivalent back plate stiffness

c2=340^2; % speed of sound squared

rho=1.205; % density of air

V=0.015; % volume of box

pref=50;

kappa=c2*rho/V;

%--------------------------------------------------------------------------

%

% Uncoupled natural frequencies of the four modes

%

%--------------------------------------------------------------------------

omt=sqrt( (Kt + kappa*At^2)/mt);

oms=sqrt( (Kb + Kt)/ms);

omb=sqrt( (Kb + kappa*Ab^2)/mb);

oma=sqrt( kappa*Aa^2/ma);

%--------------------------------------------------------------------------

%

% Coupling constants

%

%--------------------------------------------------------------------------

alphat=kappa*Aa*At;

alphab=kappa*Aa*Ab;

alphbt=kappa*Ab*At;

%--------------------------------------------------------------------------

%

% frequency, angular frequency

%

%--------------------------------------------------------------------------

f=(0:0.1:500);

om=2*pi*f;

%--------------------------------------------------------------------------

%

% Coefficients in final solution

%

%--------------------------------------------------------------------------

Dt=mt*(omt^2 - om.^2) + i*Rt*om;

Ds=ms*(oms^2 - om.^2) + i*Rs*om;

Db=mb*(omb^2 - om.^2) + i*Rb*om;

Da=ma*(oma^2 - om.^2) + i*Ra*om;

Dbar= Dt.*Ds.*Db.*Da - 2*alphbt*Kt*Kb*Da + 2*alphbt*alphab*alphat*Ds + ...

2*alphat*alphab*Kt*Kb -(alphab^2)*(Dt.*Ds - (Kt^2)) ...

-(alphat^2)*(Ds.*Db - (Kb^2)) - (alphbt^2)*Ds.*Da - (Kb^2)*Dt.*Da ...

-(Kt^2)*Db.*Da;

%--------------------------------------------------------------------------

%

% final displacements

%

%--------------------------------------------------------------------------

F=0.4; % external applied force amplitude

yt=(F./Dbar).*(Ds.*Db.*Da - (alphab^2)*Ds - (Kb^2)*Da);

ys=(F./Dbar).*(Kt*Db.*Da - Kt*(alphab^2) + alphbt*Kb*Da - alphat*alphab*Kb);

yb=(F./Dbar).*(-Kt*Kb*Da - alphbt*Da.*Ds + alphab*alphat*Ds);

ya=(F./Dbar).*(Kt*Kb*alphab + alphbt*alphab*Ds - alphat*Ds.*Db + alphat*(Kb^2) );

% The book has a minus sign in front of the first term

%--------------------------------------------------------------------------

%

% Resulting sound pressure (Equ. 2.4-11,12)

%

%--------------------------------------------------------------------------

r=1;

p=(om.*om*rho/(4*pi*r)).*(At*yt + As*ys + Ab*yb + Aa*ya);

db=20*log10(abs(p)/pref);

%--------------------------------------------------------------------------

%

% Plot results

%

%--------------------------------------------------------------------------

plot(f,db), grid on, xlabel('f (hz)'), ylabel('sound pressure level (dB)'),axis([0 500 -100 -20])