function output = Optimization_Noyrays_linear_useme(input)
% FOR OPTIMIZATION WE USE GENETIC ALGORITHM AND THE FOLLOWING DATA IN OPTIMIZATION APP
% LOWER BOUNDS [1.2 0.003 0.001 100 0.1] UPPER BOUNDS [5 0.03 0.01 500 0.2]
% INPUT 1 AIR DENSITY
% INPUT 2 THICKNESS OF THE PERFORATED PLATE
% INPUT 3 RADIUS OF THE PERFORATIONS
% INPUT 4 NUMBER OF PERFORATIONS
% INPUT 5 RADIUS 0F THE PERFORATED PLATE
rho1 = input(1); % Density of air (kg/m^3)
l = input(2); % 3mm, thickness of the plate
rp = input(3); % 1mm, radius of perforations
N = input(4); % Number of perforations
R = input(5); % 10 cm, Radius of the plate
lv = 0.1*rp;
c1 = 342; % Speed of sound (m/s)
P = N*pi*(rp^2)/(pi*R^2)
alpha = 120;
E = 7.07;
SL = 3.17*10^-1; % Laminar burning velocity
h = 0.008 ; % height of the flame monopole
A = 1.8*10^-5 ; % Flame surface area
up = 4.8; % velocity at the perforation
nn = 500;
L1 = zeros(1,nn);
for modecounter = 1:1;
switch modecounter
case 1
omegaold = complex (1000,-20);
L1 = linspace(0.1,0.8,nn);
case 2
omegaold = complex (2200,-60);
L1 = zeros(1,nn);
L1 = linspace(0.1,0.8,nn);
case 3
omegaold = complex (3500,130);
L1 = zeros(1,nn);
L1 = linspace(0.1,0.8,nn);
end
for nL1 = nn:-1:1;
Norma_amp = 0;
n1 = -1.185*(10^-1)*Norma_amp + 1.99;
n2 = 1.185*(10^-1)*Norma_amp - 0.99;
tau1 = 2.28*(10^-3)*(Norma_amp^2) - 4.74*(10^-5)*Norma_amp + 8.7*(10^-4);
tau2 = 4.35*(10^-3)*(Norma_amp^2) + 6.79*(10^-5)*Norma_amp + 1.12*(10^-3);
sigma1 = 5.86*(10^-4)*Norma_amp + 2.36*(10^-4);
sigma2 = -4.63*(10^-4)*Norma_amp + 9.3*(10^-4);
w0(nL1) = w0optimization_secant_method_Saad_Jan(L1(nL1),omegaold,c1,rho1,l,rp,lv,N,R,P); % This function is given at the end of this code
omegaold = w0(nL1);
w0r(nL1) = real(w0(nL1));
w0i(nL1) = imag(w0(nL1));
FDF(nL1) = n1*exp(-(w0r(nL1))^2 *0.5*sigma1^2)*exp(i*(w0r(nL1))*tau1) + n2*exp(-(w0r(nL1))^2 *0.5*sigma2^2)*exp(i*(w0r(nL1))*tau2);
Gainvect(nL1) = abs(FDF(nL1));
Gain = Gainvect(nL1);
Phase(nL1) = angle(FDF(nL1));
taub(nL1) = L1(nL1)/c1; % characteristic acoustic time of the burner
taup(nL1) = l/(c1*P); % characteristic acoustic time delay of the plate
tauco(nL1) = alpha*(E-1)*SL*A*Gain/(up*P*4*pi*h*c1); % acoustic flame coupling timescale
taucv = 1.26*10^-3;
w0star(nL1) = w0(nL1)*taub(nL1);
epsilon(nL1) = tauco(nL1)/taub(nL1);
w1star(nL1) = -epsilon(nL1)*w0star(nL1)*exp(i*w0star(nL1)*taucv/taub(nL1))/(1 + (1 + w0star(nL1)^2*taup(nL1)/taub(nL1))*taup(nL1)/taub(nL1));
wr(nL1) = w0r(nL1)*(1-epsilon(nL1)*cos(w0r(nL1)*taucv)/((1 + (w0r(nL1))^2*taup(nL1)*taub(nL1))*(1+taup(nL1)/taub(nL1))));
wi(nL1) = -epsilon(nL1)*w0r(nL1)*sin(w0r(nL1)*taucv)/(1+(taup(nL1)/taub(nL1))*(1 + (w0r(nL1))^2*taup(nL1)*taub(nL1)));
w(nL1) = wr(nL1) + i*wi(nL1);
f(nL1) = wr(nL1)/(2*pi);
end
switch modecounter
case 1
wimax1 = max(wi);
output = wimax1;
case 2
wimax2 = max(wi);
output = wimax2;
case 3
wimax3 = max(wi);
output = wimax3;
end
end
% The other function w0optimization_secant_method_Saad_Jan(L1(nL1),omegaold,c1,rho1,l,rp,lv,N,R,P) which finds the value of w0 in the code above is given below
function omega = w0optimization_secant_method_Saad_Jan(L1,omegaold,c1,rho1,l,rp,lv,N,R,P)
error = 1;
counter = 0;
omega = omegaold;
omegaold = 0.95*omega;
while error > 1e-9 %1e-9
counter = counter+1;
F = omega*rho1*l*(1 + lv*rp*(1+1i)) - rho1*c1*P*cot(omega*L1/c1); % Eq. defining the case for no combustion
Fold = omegaold*rho1*l*(1 + lv*rp*(1+1i)) - rho1*c1*P*cot(omegaold*L1/c1);
omeganew = omega - F*(omega-omegaold)/(F-Fold);
error = abs((omeganew-omega)/omeganew);
omegaold = omega;
omega = omeganew;
end
The error pic of the optimization tool box is attached as below
