کد متلب مدار شانت پیزوالکتریک

به درخواست یکی از بازدید کننده ها قسمتی از کد متلب نوشته شده برای مدار شانت پیزوالکتریک بر روی تیر رو گذاشتم. امید وارم به کمکتون بیاد.

کد نوشته شده در متلب

% PARALLEL R-L SHUNT OF HINGED-HINGED BEAM

% WRITTEN BY A.Keshavarz

% ar.keshavarz.aut.ac.ir

CLEAR

B=0.025; % BEAM & PZT WIDTH

L=0.5; % BEAM LENGTH

X1=0.095; % DISTANCE OF THE LEFT EDGE OF PZT

X2=0.155; % DISTANCE OF THE RIGHT EDGE OF PZT, (X2-X1) IS THE LENGTH OF PZT

XD=0.5*L; % POSITION OF FORCE APPLIED

EB=70E9; % YOUNG'S MODULUS OF BEAM

DB=2700; % DENSITY OF BEAM

HB=3E-3; % THICKNESS OF BEAM

EP=70E9; % YOUNG'S MODULUS OF PZT

DP=7600; % DENSITY OF PZT

HP=0.25E-3; % THICKNESS OF PZT

H31=7.65E8; % PZT CONSTANT

AB=B*HB; % AREA OF BEAM

AP=B*HP; % AREA OF PZT

ZB=HB/2; % CENTER LOCATION OF BEAM

ZP=HB+HP/2; % CENTER LOCATION OF PZT

ZN=(AB*EB*ZB+AP*EP*ZP)/(AB*EB+AP*EP); % NATURAL AXIS OF COMPOSITE BEAM

IB=B*HB^3/12+AB*(ZN-ZB)^2; % MOMENT OF INERTIA OF BEAM

IP=B*HP^3/12+AP*(ZN-ZP)^2; % MOMENT OF INERTIA OF PZT

N=3; % NO. OF FREQUENCIES

% INITIALIZE MATRIX

K=ZEROS(N);

M=ZEROS(N);

C=ZEROS(N);

BL=ZEROS(N,1);

BR=ZEROS(N,1);

FD=ZEROS(N,1);

CP=82E-9; % CAPACITANCE OF PZT

%L0=[400 4.788 1]; % SHUNT INDUCTORS

%R0=[3E5 1E4 1]; % SHUNT RESISTORS

L0=[500 24.64 4.846]; % SHUNT INDUCTORS

R0=[3.4E5 1E5 7E4]; % SHUNT RESISTORS

% COEFFICIENT OF BL

%BL=H31*(HB+0.5*HP)*(L0/R0)*(PI/L)*(COS(PI*X2/L)-COS(PI*X1/L));

% COEFFICIENT OF BR

%BR=H31*(HB+0.5*HP)*(PI/L)*(COS(PI*X2/L)-COS(PI*X1/L));

% COEFFICIENT OF BL

%BL=H31*(HP/(X2-X1))*(HB+0.5*HP)*(L0/R0)*(PI/L)*...

%(COS(PI*X2/L)-COS(PI*X1/L));

% COEFFICIENT OF BR

%BR=H31*(HP/(X2-X1))*(HB+0.5*HP)*(PI/L)*(COS(PI*X2/L)-COS(PI*X1/L));

CR=INV(DIAG(R0*CP)); % COEFFICIENT OF CR

CL=INV(DIAG(L0*CP)); % COEFFICIENT OF CL

FOR R=1:N;

Q1(R)=-(H31*HP*(HB+0.5*HP)*(R*PI/L)*ABS(COS(R*PI*X2/L)-...

COS(R*PI*X1/L)))/((X2-X1));

FOR S=1:N;

BL(S)=H31*(HP/(X2-X1))*(HB+0.5*HP)*(L0(S)/R0(S))*...

(S*PI/L)*(COS(PI*X2/L)-COS(PI*X1/L)); % COEFFICIENT OF BL

BR(S)=H31*(HP/(X2-X1))*(HB+0.5*HP)*(S*PI/L)*...

(COS(PI*X2/L)-COS(PI*X1/L)); % COEFFICIENT OF BR

IF R == S

K(R,S)=(PI*R/L)^4*(EB*IB*L/2+EP*IP*(X2-X1)/2+...

EP*IP*L/(4*PI*R)*(SIN(2*PI*R*X1/L)-SIN(2*PI*R*X2/L)));

M(R,S)=(DB*AB*L/2+DP*AP*(X2-X1)/2+DP*AP*L/...

(4*PI*R)*(SIN(2*PI*R*X1/L)-SIN(2*PI*R*X2/L)));

ELSE

K(R,S)=EP*IP*L/PI*(PI^2*R*S/L^2)^2*((R*SIN(S*PI*X2/L)...

*COS(R*PI*X2/L))/(S^2-R^2)+(S*SIN(R*PI*X2/L)*...

COS(S*PI*X2/L))/(R^2-S^2)-((R*SIN(S*PI*X1/L)*...

COS(R*PI*X1/L))/(S^2-R^2)+(S*SIN(R*PI*X1/L)*...

COS(S*PI*X1/L))/(R^2-S^2)));

M(R,S)=DP*AP*L/PI*((R*SIN(S*PI*X2/L)*COS(R*PI*X2/L))/...

(S^2-R^2)+(S*SIN(R*PI*X2/L)*COS(S*PI*X2/L))/(R^2-S^2)-...

((R*SIN(S*PI*X1/L)*COS(R*PI*X1/L))/(S^2-R^2)+...

(S*SIN(R*PI*X1/L)*COS(S*PI*X1/L))/(R^2-S^2)));

END;

END;

%FD(R)=1/100*SIN(R*PI*XD/L); % DUE TO DISCRETE FORCE WITH MAGNITUDE 1/100

END;

FD = (1/100)*[1; 1; 1]; % APPLIED FORCE

CD=0.2*M+1E-7*K; % INTERNAL DAMPING COEFFICIENT

Q=Q1*INV(DIAG(L0)); % Q1/L0

Q=DIAG(Q); % MAKE Q AS A DIAGONAL MATRIX

BL=DIAG(BL);

BR=DIAG(BR);

% STATE-SPACE MODEL

AK=-INV(M)*K; % -K/M

AC=-INV(M)*CD; % -C/M

% A MATRIX WITHOUT SHUNT

A1=[ZEROS(N) EYE(N) ZEROS(N) ZEROS(N); AK AC ZEROS(N) ZEROS(N); ...

ZEROS(N) ZEROS(N) ZEROS(N) EYE(N); -Q ZEROS(N) -CL -CR];

% A MATRIX WITH SHUNT

A2=[ZEROS(N) EYE(N) ZEROS(N) ZEROS(N); AK AC -INV(M)*BR -INV(M)*BL; ...

ZEROS(N) ZEROS(N) ZEROS(N) EYE(N); -Q ZEROS(N) -CL -CR];

B=[ZEROS(N,1); INV(M)*FD; ZEROS(N,1); ZEROS(N,1)];

FOR R=1:N;

C1(1,R)=SIN(R*PI/4); % DISPLACEMENT W AT MIDPOINT (X=L/2)

END;

C=[C1 ZEROS(1,3*N)];

D=[0];

% IMPULSE RESPONSE

T=LINSPACE(0,1,1E5);

IU=1;

[Y,X,T]=IMPULSE(A1,B,C,D,IU,T); % UNCONTROLLED RESPONSE

[YC,X,T]=IMPULSE(A2,B,C,D,IU,T); % CONTROLLED RESPONSE

FIGURE(1);

PLOT(T,Y,'B--',T,YC,'R') % UNIT (MM)

LEGEND('WITHOUT PZT SHUNT','WITH PZT SHUNT');

TITLE('IMPULSE RESPONSE OF BEAM TRANSVERSE DISPLACEMENT AT XM=L/4', 'FONTSIZE', 18, 'COLOR', 'R')

XLABEL('TIME (SEC)', 'FONTSIZE', 16, 'COLOR', 'R')

YLABEL('DISPLACEMENT (M)', 'FONTSIZE', 16, 'COLOR', 'R')

AXIS([0 1 -1E-3 1E-3])

GRID;

ZOOM ON;

FIGURE(2);

PLOT(20*LOG10(ABS(FFT(Y))),'B--');

HOLD ON;

PLOT(20*LOG10(ABS(FFT(YC))), 'R');

LEGEND('WITHOUT PZT SHUNT','WITH PZT SHUNT');

TITLE('FREQUENCY RESPONSE OF BEAM TRANSVERSE DISPLACEMENT AT XM=L/4','FONTSIZE', 20, 'COLOR', 'R');

XLABEL('FREQUENCY (HZ)', 'FONTSIZE', 16, 'COLOR', 'R');

YLABEL('MAGNITUDE (DB)', 'FONTSIZE', 16, 'COLOR', 'R');

AXIS([0 1000 -80 40]);

GRID;

ZOOM ON;

HOLD OFF;