Контекст:
Проблема, которую я решаю, состоит в том, чтобы найти оптимальное размещение и размер пьезоэлектрического c патча на балке, чтобы модальная сила была максимальной. В этом коде я вычисляю модальные формы с помощью метода Ритца, а затем применяю уравнение, чтобы получить модальную силу, а затем суммировать различные режимы и делить на количество режимов, которые я добавляю. Цель ограничения - минимизировать дрейф собственных частот (то есть представить очень тонкий и короткий луч с большим участком пьезоэлектрического c - не очень хорошо). Я реализую это путем ограничения всех модифицированных собственных частот в пределах n-мерного шара с радиусом c* величина (омега), где c - число от 0 до 1.
Код :
Point Density
M = 500;
H = linspace(0.075e-3,1e-2,M);
E = 200*10^9;
rho = 7700;
xi_opt = zeros(2,M);
for j = 1:M
[xi_opt(:,j)] = Optimizer(H(j),E,rho);
end
function [xi] = Optimizer(h,E,rho)
% Constraint value
c = 0.3;
L = 1;
xi_0 = [0.2;0.01];
A = [1,2];
b = [L];
Aeq = [];
beq = [];
lb = eps.*ones(2,1);
ub = Inf.*ones(2,1);
options = optimoptions('patternsearch','display','iter',...
'ConstraintTolerance',1e-8,'FunctionTolerance',1e-3...
,'MaxFunctionEvaluations',10000,'MeshTolerance',1e-3,'StepTolerance'...
,1e-3);
[xi] = patternsearch(@(xi) f(xi,h,E,rho),xi_0,A,b,Aeq,beq,lb,ub,@(xi)...
nonlcon(xi,c,h,E,rho),options);
end
function [ciq,ceq] = nonlcon(xi,c,h,E,rho)
[~,g] = f(xi,h,E,rho);
ciq = g-c;
ceq = [];
end
function [metric1,metric2] = f(xi,h,E,rho)
%% Material Properties
% Beam properties
L=1; % beam length (m)
b=0.1; % beam width (m)
%{
h=.24e-3; % thickness (m)
E=200*10^9; % Young's modulus
rho=7700; % density
%}
% Properties of the Piezoelectric Film
E_p=30.336e9; %Young's modulus (Pa)
d33=460e-12; %piezo constant, elongating type (C/N)
h_p=0.30e-3; %thickness of MFC, (m)
l_p=0.45e-3; %distance between poles (m)
rho_p=5440; %density (kg/m^3)
%% Modal settings
N=20; % number of Ritz vectors to use (with 13, first five non-zero
% % modes are well converged)
epw=5000; % number of elements per half-wave to discretize beam into
% Voltage amplitude
V=400; %Amplitude (+/-400 is about the maximum safe limit)
%% Piezoelectric Calcs
% Distance from neutral axis to geometric mid-plane
y0=(E_p*h_p*(h+h_p))/(2*(E*h+E_p*h_p));
% Electromechanical coupling coefficient
chi=(d33*E_p*b*h_p/l_p)*((.5*(h+h_p)-y0));
% No adjustments for patch stiffness or mass
[wn, ~, ~]=DDSL_Beam_Modes(L,b,h,E,rho,N,epw);
% Adjustments for patch stiffness or mass
[wnp, Phip, nodemap]=DDSL_Beam_Modes_p(L,b,h,h_p,E,E_p,rho,...
rho_p,N,epw,xi);
%% Calcs on modal outputs
% Initialize secondary counters for average and distance
g = 0;
h = 0;
metric1 = 0;
metric2 = 0;
% Initialize modal force vector
Q = zeros(N,1);
% differential distance
dx = nodemap(2)-nodemap(1);
% Node points at beginning and end of piezoelectric element
k = dsearchn(nodemap',[xi(2)*L; L*(xi(1)+xi(2))]);
% Compute modal force (See Erturk and Inman Book Eq. 3.31)
for n = 1:N
dphi=gradient(Phip(:,n),dx); % Use piezo-adjusted modes
Q(n)=chi*V*(dphi(k(2))-dphi(k(1)));
if (n>2)&&(n<8)
metric1 = metric1+abs(Q(n));
metric2 = metric2+((wn(n)-wnp(n)))^2;
g = g+1;
h = h+(wn(n))^2;
end
end
% Normalize metrics
metric1 = -(metric1/g);
metric2 = sqrt((metric2/h));
end
Я также добавлю код для DDSL_Beam_Modes_p. Код для DDSL_Beam_Modes очень похож, но контуры, определяющие M и K, немного отличаются, поскольку выражение для безразмерного центра тяжести изменено из-за присутствия элемента piezoelectri c. Я могу предоставить применимые теоретические основы этого и метода Ритца, если необходимо.
function [wnp, Phip, nodemap]=DDSL_Beam_Modes_p(L,b,h,hp,E,Ep,rho,...
rhop,N,epw,xi)
% Finds the natural frequencies and modes of a FF beam using the Ritz
% approach.
% Includes effects of piezo stiffness and mass.
%% Define and Initialize Variables
xg=L/2; %get center point of plate (for FF without piezo)
xig = (0.5*h*rho+(xi(2)+0.5*xi(1))*(xi(1)*hp*rhop))/(h*rho+hp*xi(1)*rhop);
nx=epw*N; %discretize into epw times the highest mode number
YI=E*((b*(h^3))/12); %Flexural rigidity of substrate
cx=linspace(0,L,nx);
nodemap=cx;
% Piezo calculations - calculate height and moment of inertia in segment of
% beam with piezoelectric patch
y0=(Ep*hp*(h+hp))/(2*(E*h+Ep*hp));
YIc=b*(E*(h^3/12+h*y0^2)+Ep*(hp^3/12+hp*(.5*(h+hp)-y0)^2));
% Admissible function for FF beam
psi = @(x,j) (((x-xg)./L).^(j-1));
%% Calculate M and K via Ritz Method
M=zeros(N,N);
K=zeros(N,N);
for j=1:N
for i=1:N
M(i,j)=rho*h*b*((L/(i+j-1))*((-xig+xi(2))^(i+j-1)-(-xig)^...
(i+j-1)))+(rho*h+rhop*hp)*b*((L/(i+j-1))*((-xig+xi(1)+...
xi(2))^(i+j-1)-(-xig+xi(2))^(i+j-1)))+rho*h*b*((L/(i+j-1))*...
((xig)^(i+j-1)-(-xig+xi(2)+xi(1))^(i+j-1)));
end
end
for j=3:N
for i=3:N
K(i,j)=((((j-1)*(j-2)*(i-1)*(i-2))/(L^3))*((1)/(i+j-5)))*...
(YI*((-xig+xi(2))^(i+j-5)-(-xig)^(i+j-5))+YIc*(((-xig+xi(1)+...
xi(2))^(i+j-5)-(-xig+xi(2))^(i+j-5)))+...
YI*((xig)^(i+j-5)-(-xig+xi(2)+xi(1))^(i+j-5)));
end
end
%% Apply Classical Modal Theory
% Find eigenfrequencies
[phi,lambda] = eig(K,M);
[wnp,ind_sort] = sort(sqrt(diag(lambda)));
phi_sort = zeros(N,N);
Phim = zeros(N,N);
for j=1:N
phi_sort(:,j) = phi(:,ind_sort(j));
%Mass normalize mode shapes
Phim(:,j) = phi_sort(:,j)./sqrt(phi_sort(:,j)'*M*phi_sort(:,j));
end
Psi = zeros(nx,N);
for j = 1:N
Psi(:,j) = psi(cx,j);
end
Phip = Psi*Phim; %Transform into physical coordinates
end
Графики и рисунки:
Для полноты картины я включу графики f и g для конкретной толщины балки. Надеюсь, это будет полезно.
Графики f и g
Вопрос:
Когда N (количество функций Ритца) равно 24, код работает невероятно медленно, но работает успешно, но когда я увеличиваю его до N = 25, появляется следующая ошибка:
>> DDSL_Piezo_Optimizer
Max
Iter Func-count f(x) Constraint MeshSize Method
0 1 -2.90857e+06 1.371 1
1 15 -2.90857e+06 1.371 0.001 Increase penalty
Operands to the || and && operators must be convertible to logical scalar values.
Error in psAugConverged (line 94)
if ~isempty(infMessage) && strmatch('optimlib:optimfcnchk',infMessage)
Error in pfmincon (line 66)
[X,FVAL,maxConstr,optimState] = psAugConverged(Iterate,X,optimState,psAugParam,options);
Error in patternsearch (line 280)
[X,FVAL,EXITFLAG,OUTPUT] = pfmincon(FUN,X0,optimState,Iterate, ...
Error in DDSL_Piezo_Optimizer>Optimizer (line 32)
[xi] = patternsearch(@(xi) f(xi,h,E,rho),xi_0,A,b,Aeq,beq,lb,ub,@(xi)...
Error in DDSL_Piezo_Optimizer (line 11)
[xi_opt(:,j)] = Optimizer(H(j),E,rho);
>>
Я понимаю, что вызывает эту ошибку - пытаюсь сравнить массивы с && и || вместо & и |, но я не понимаю, почему эта ошибка не существует, когда N <25. </p>