%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Nonlinear Simulink Model for a MEMS Microphone - Companion Script %
% Max Lacoma, 2024 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Details
%
% This is a companion script for two Simulink models (nonlinear_mems_mic
% and nonlinear_ported_mems_mic).
%
% The first model explores the behavior of a MEMS microphone when certain
% nonlinearities have been included in the model. The first nonlinearity is
% the acoustic resistance due to the viscous forces in the gap between the
% membrane and back plate. This resistance is dependent on the diaphragm
% displacement, which is dynamically computed in this model. Secondly, the
% force and capacitance of the electrostatic gap are dependent on both the
% diaphragm displacement and the voltage across the gap. Both of these are
% used to compute the variable capacitance and force for the model.
%
% Various scopes have been added to the model to explore the
% nonlinearities. Additionally, running this script will show the
% linearized approximation of the frequency response at varying bias
% voltages. The linear part of the model and all parameters were found in
% "Structure-based equivalent circuit modeling of a capacitive-type MEMS
% microphone" by Lee et. al. The model was then expanded upon by Max Lacoma
% to include the nonlinearities.
%
% The second model looks at the effect of a microphone port on the
% frequency response. Specifically, it models the outer case, gasket,
% acoustic mesh, and PCB hole for a porting scheme with increasing
% diameter. Realistic geometry for this porting scheme was found in "Finite
% Element Modeling of MEMS Microphone Ports" by Jerad Lewis. Running this
% script will show the linearized frequency responses of the non-ported and
% ported geometries.
%% Set Parameters for Models
clear; close all; clc
%Frequency for computing the frequency dependent radiation resistance at
f = 20000;
w = 2*pi*f;
%Geometric parameters
r = 300e-6; %Diaphragm radius [m]
N = 1668; %Number of backplate holes
tb = 2.5e-6; %Backplate thickness [m]
rb = 4e-6; %Radius of backplate holes [m]
Ap = (pi*rb^2*N)/(pi*r^2); %Fraction of backplate covered in holes
Vbc = 6e-9; %Back volume [m^3]
Deff = 0.6e-3; %Effective diameter of diaphragm [m]
Aeff = Deff^2*pi/4; %Effective area of diaphragm [m^2]
t = 0.9e-6; %Diaphragm thickness [m]
d0 = 2.45e-6; %Diaphragm displacement at zero bias [m]
%Material parameters
rho = 1.21; %Air density [kg/m^3]
nu = 1.9841e-05; %Dynamic viscocity of air [kg/(m*s)]
c = 343; %Speed of sound in air [m/s]
rhod = 3809; %Equivalent density of diaphragm [kg/m^3]
%Acoustic domain parameters
%Note that Rg is nonlinear - needs d as an input - gap separation
Mr = 8*rho/(3*pi^2*r); %Acoustic radiation inductance [kg/m^4]
Rr = rho*w^2/(2*pi*c); %Acoustic radiation resistance [kg/(s*m^4)]
Rh = 8*nu*tb/(N*pi*rb^4); %Acoustic resistance of backplate holes [kg/(s*m^4)]
Rgtimesd3 = 12*nu/(N*pi)* (Ap/2-Ap^2/8-log(Ap)/4-3/8); %Acoustic resistance of backplate gap [kg/(s*m^4)]
Cbc = Vbc/(rho*c^2); %Acoustic compliance of rear volume [m^3/Pa]
%Mechanical domain parameters
Md = 78.56e-11; %Effective mass of the diaphragm [kg]
Cd = 3.43e-3; %Effective compliance of diaphragm [m/N]
%Electrical domain parameters
Cp = 0.2e-12; %Parasitic capacitance of microphone [F]
e0 = 8.85e-12; %Vacuum permittivity [F/m]
%% Plot Frequency Responses @ Varying Bias Voltages - Non-Ported
%Set bias voltage
Vbias_vec = [3.3 5 10.4]; %[V]
%Define input and output points for linearization
io(1) = linio('nonlinear_mems_mic/1 kHz 1 Pa Input',1,'input');
io(2) = linio('nonlinear_mems_mic/Vout',1,'output');
%Create a figure
figure; hold on
%Loop through each bias voltage
for ii = 1:numel(Vbias_vec)
%Set bias voltage in model
Vbias = Vbias_vec(ii);
%Linearize the model at 1 second
sys = linearize('nonlinear_mems_mic',io,1);
%Compute the frequency response
[mag,phase,wout] = bode(sys);
%Parse system
mag2 = zeros(1,numel(mag));
phase2 = mag2;
for ii = 1:numel(mag)
mag2(ii) = mag(1,1,ii);
phase2(ii) = phase(1,1,ii);
end
clear mag phase
%Plot in dBV/Pa
semilogx(wout/(2*pi),20*log10(mag2))
end
%Set plot settings
xlabel('Frequency [Hz]')
ylabel('Magnitude [dBV/Pa, re: 1 V]')
leg = legend('V_b_i_a_s = 3.3 V','V_b_i_a_s = 5.0 V','V_b_i_a_s = 10.4 V');
title('Linearized MEMS Mic Frequency Response')
axis tight
set(gca,'XScale','log')
leg.Location = 'northwest';
%% Port Parameters
%Geometry of ports
a_pcb = 0.8e-3; %PCB hole radius [m]
L_pcb = 0.5e-3; %PCB hole length [m]
a_gasket = 1e-3; %Gasket hole radius [m]
L_gasket = 1.6e-3; %Gasket hole length [m]
a_case = 1e-3; %Case hole radius [m]
L_case = 4e-3; %Case hole length [m]
%Acoustic masses and compliances for the ports
[m_pcb C_pcb ~] = acs_mass_comp(a_pcb,L_pcb,rho,c);
[m_gasket C_gasket ~] = acs_mass_comp(a_gasket,L_gasket,rho,c);
[m_case C_case A_case] = acs_mass_comp(a_case,L_case,rho,c);
%Acoustic mesh
Rmesh = 40; %[Rayl]
Rmesh = Rmesh/a_gasket; %[Rayl/m^2]
%% Compare Ported vs Non-Ported Frequency Responses
%Define input and output points for linearization
io(1) = linio('nonlinear_ported_mems_mic/1 kHz 1 Pa Input',1,'input');
io(2) = linio('nonlinear_ported_mems_mic/Vout',1,'output');
%Create a figure - plot 10.4V response of nonported mic
figure; hold on
semilogx(wout/(2*pi),20*log10(mag2))
%Linearize the model at 1 second
sys = linearize('nonlinear_ported_mems_mic',io,1);
%Compute the frequency response
[mag,phase,wout] = bode(sys);
%Parse system
mag2 = zeros(1,numel(mag));
phase2 = mag2;
for ii = 1:numel(mag)
mag2(ii) = mag(1,1,ii);
phase2(ii) = phase(1,1,ii);
end
clear mag phase
%Plot in dBV/Pa
semilogx(wout/(2*pi),20*log10(mag2))
%Set plot settings
xlabel('Frequency [Hz]')
ylabel('Magnitude [dBV/Pa, re: 1 V]')
leg = legend('Non-Ported','Ported');
title('Linearized MEMS Mic Frequency Response')
axis tight
xlim([0 200000])
set(gca,'XScale','log')
leg.Location = 'northwest';
%% Functions
%Function to compute the acoustic mass and compliance of the ports
function [m C A] = acs_mass_comp(a,L,rho,c)
%Compute area and volume
A = pi*a^2;
V = A*L;
%Compute mass and compliance
m = rho*L/A;
C = V/(rho*c^2);
end