📄 act8pz.m
字号:
echo off
% act8pz.m Matlab model of two-stage actuator/suspension system,
% balanced reduction.
clear all;
hold off;
clf;
% load the Block Lanczos .mat file actrl_eig.mat, containing evr - the modal matrix, freqvec -
% the frequency vector and node_numbers - the vector of node numbers for the modal
% matrix
% the output for the ANSYS run is the following dof's
% dof node dir where
%
% 1 22 ux - radial, top head gap
% 2 10022 ux - radial, bottom head gap
% 3 21538 ux - radial, bottom arm piezo, hub end
% 4 21546 ux - radial, bottom arm piezo, head end
% 5 21576 ux - radial, bottom arm piezo, hub end
% 6 21584 ux - radial, bottom arm piezo, head end
% 7 21617 ux - radial, bottom arm piezo, hub end
% 8 21625 ux - radial, bottom arm piezo, head end
% 9 22538 ux - radial, bottom arm piezo, hub end
% 10 22546 ux - radial, bottom arm piezo, head end
% 11 22576 ux - radial, bottom arm piezo, hub end
% 12 22584 ux - radial, bottom arm piezo, head end
% 13 22617 ux - radial, bottom arm piezo, hub end
% 14 22625 ux - radial, bottom arm piezo, head end
% 15 24061 ux - radial, bottom arm piezo, coil
% 16 24066 ux - radial, bottom arm piezo, coil
% 17 24082 ux - radial, bottom arm piezo, coil
% 18 24087 ux - radial, bottom arm piezo, coil
% 19 24538 ux - radial, top arm piezo, hub end
% 20 24546 ux - radial, top arm piezo, head end
% 21 24576 ux - radial, top arm piezo, hub end
% 22 24584 ux - radial, top arm piezo, head end
% 23 24617 ux - radial, top arm piezo, hub end
% 24 24625 ux - radial, top arm piezo, head end
% 25 25538 ux - radial, top arm piezo, hub end
% 26 25546 ux - radial, top arm piezo, head end
% 27 25576 ux - radial, top arm piezo, hub end
% 28 25584 ux - radial, top arm piezo, head end
% 29 25617 ux - radial, top arm piezo, hub end
% 30 25625 ux - radial, top arm piezo, head end
% 31 22 uy - circumferential, top head gap
% 32 10022 uy - circumferential, bottom head gap
% 33 21538 uy - circumferential, bottom arm piezo, hub end
% 34 21546 uy - circumferential, bottom arm piezo, head end
% 35 21576 uy - circumferential, bottom arm piezo, hub end
% 36 21584 uy - circumferential, bottom arm piezo, head end
% 37 21617 uy - circumferential, bottom arm piezo, hub end
% 38 21625 uy - circumferential, bottom arm piezo, head end
% 39 22538 uy - circumferential, bottom arm piezo, hub end
% 40 22546 uy - circumferential, bottom arm piezo, head end
% 41 22576 uy - circumferential, bottom arm piezo, hub end
% 42 22584 uy - circumferential, bottom arm piezo, head end
% 43 22617 uy - circumferential, bottom arm piezo, hub end
% 44 22625 uy - circumferential, bottom arm piezo, head end
% 45 24061 uy - circumferential, bottom arm piezo, coil
% 46 24066 uy - circumferential, bottom arm piezo, coil
% 47 24082 uy - circumferential, bottom arm piezo, coil
% 48 24087 uy - circumferential, bottom arm piezo, coil
% 49 24538 uy - circumferential, top arm piezo, hub end
% 50 24546 uy - circumferential, top arm piezo, head end
% 51 24576 uy - circumferential, top arm piezo, hub end
% 52 24584 uy - circumferential, top arm piezo, head end
% 53 24617 uy - circumferential, top arm piezo, hub end
% 54 24625 uy - circumferential, top arm piezo, head end
% 55 25538 uy - circumferential, top arm piezo, hub end
% 56 25546 uy - circumferential, top arm piezo, head end
% 57 25576 uy - circumferential, top arm piezo, hub end
% 58 25584 uy - circumferential, top arm piezo, head end
% 59 25617 uy - circumferential, top arm piezo, hub end
% 60 25625 uy - circumferential, top arm piezo, head end
load actrlpz_eig;
[numdof,num_modes_total] = size(evr);
freqvec(1) = 0; % set rigid body mode to zero frequency
xn = evr;
% define radial and circumferential forces applied at four coil force nodes
% "x" is radial, "y" is circumferential, total force is unity
fcoil = 0.25;
n24061fx = fcoil*sin(9.1148*pi/180);
n24061fy = fcoil*cos(9.1148*pi/180);
n24066fx = fcoil*sin(15.1657*pi/180);
n24066fy = fcoil*cos(15.1657*pi/180);
n24082fx = -fcoil*sin(15.1657*pi/180);
n24082fy = fcoil*cos(15.1657*pi/180);
n24087fx = -fcoil*sin(9.1148*pi/180);
n24087fy = fcoil*cos(9.1148*pi/180);
% define radial and circumferential forces applied at ends of piezo element
% "x" is radial, "y" is circumferential, total force is unity
fpz = 0.2/6; % six nodes at each end of the piezo
% bottom arm radial force
n21538fx = fpz*cos(20.4549*pi/180);
n21546fx = -fpz*cos(13.5298*pi/180);
n21576fx = fpz*cos(20.4549*pi/180);
n21584fx = -fpz*cos(13.5298*pi/180);
n21617fx = fpz*cos(20.4549*pi/180);
n21625fx = -fpz*cos(13.5298*pi/180);
n22538fx = fpz*cos(20.4549*pi/180);
n22546fx = -fpz*cos(13.5298*pi/180);
n22576fx = fpz*cos(20.4549*pi/180);
n22584fx = -fpz*cos(13.5298*pi/180);
n22617fx = fpz*cos(20.4549*pi/180);
n22625fx = -fpz*cos(13.5298*pi/180);
% top arm radial force
n24538fx = -fpz*cos(20.4549*pi/180);
n24546fx = fpz*cos(13.5298*pi/180);
n24576fx = -fpz*cos(20.4549*pi/180);
n24584fx = fpz*cos(13.5298*pi/180);
n24617fx = -fpz*cos(20.4549*pi/180);
n24625fx = fpz*cos(13.5298*pi/180);
n25538fx = -fpz*cos(20.4549*pi/180);
n25546fx = fpz*cos(13.5298*pi/180);
n25576fx = -fpz*cos(20.4549*pi/180);
n25584fx = fpz*cos(13.5298*pi/180);
n25617fx = -fpz*cos(20.4549*pi/180);
n25625fx = fpz*cos(13.5298*pi/180);
% bottom arm circumferential force
n21538fy = -fpz*sin(20.4549*pi/180);
n21546fy = fpz*sin(13.5298*pi/180);
n21576fy = -fpz*sin(20.4549*pi/180);
n21584fy = fpz*sin(13.5298*pi/180);
n21617fy = -fpz*sin(20.4549*pi/180);
n21625fy = fpz*sin(13.5298*pi/180);
n22538fy = -fpz*sin(20.4549*pi/180);
n22546fy = fpz*sin(13.5298*pi/180);
n22576fy = -fpz*sin(20.4549*pi/180);
n22584fy = fpz*sin(13.5298*pi/180);
n22617fy = -fpz*sin(20.4549*pi/180);
n22625fy = fpz*sin(13.5298*pi/180);
% top arm circumferential force
n24538fy = fpz*sin(20.4549*pi/180);
n24546fy = -fpz*sin(13.5298*pi/180);
n24576fy = fpz*sin(20.4549*pi/180);
n24584fy = -fpz*sin(13.5298*pi/180);
n24617fy = fpz*sin(20.4549*pi/180);
n24625fy = -fpz*sin(13.5298*pi/180);
n25538fy = fpz*sin(20.4549*pi/180);
n25546fy = -fpz*sin(13.5298*pi/180);
n25576fy = fpz*sin(20.4549*pi/180);
n25584fy = -fpz*sin(13.5298*pi/180);
n25617fy = fpz*sin(20.4549*pi/180);
n25625fy = -fpz*sin(13.5298*pi/180);
% two-input system
% first input is coil force
% second input is excitation of both piezo elements with opposite polarity
% f_coil is the vector of forces applied to coil
f_coil = [zeros(14,1)
n24061fx
n24066fx
n24082fx
n24087fx
zeros(26,1)
n24061fy
n24066fy
n24082fy
n24087fy
zeros(12,1)];
% f_piezo is vector of forces applied to piezo ends
f_piezo = [ 0
0
n21538fx % bottom arm radial force
n21546fx
n21576fx
n21584fx
n21617fx
n21625fx
n22538fx
n22546fx
n22576fx
n22584fx
n22617fx
n22625fx
0
0
0
0
n24538fx % top arm radial force
n24546fx
n24576fx
n24584fx
n24617fx
n24625fx
n25538fx
n25546fx
n25576fx
n25584fx
n25617fx
n25625fx
0
0
n21538fy % bottom arm circumferential force
n21546fy
n21576fy
n21584fy
n21617fy
n21625fy
n22538fy
n22546fy
n22576fy
n22584fy
n22617fy
n22625fy
0
0
0
0
n24538fy % top arm circumferential force
n24546fy
n24576fy
n24584fy
n24617fy
n24625fy
n25538fy
n25546fy
n25576fy
n25584fy
n25617fy
n25625fy ];
% define composite forcing function, force applied to each node times eigenvector value
% for that node
force_coil = f_coil'*xn;
force_piezo = f_piezo'*xn;
% prompt for uniform or variable zeta
zeta_type = input('enter "1" to read in damping vector (zetainpz.m) or "enter" for uniform damping ... ');
if (isempty(zeta_type))
zeta_type = 0;
zeta_uniform = input('enter value for uniform damping, .005 is 0.5% of critical (default) ... ');
if (isempty(zeta_uniform))
zeta_uniform = 0.005;
end
zeta_unsort = zeta_uniform*ones(num_modes_total,1);
else
zetainpz; % read in zeta_unsort damping vector from zetainpz.m file
end
if length(zeta_unsort) ~= num_modes_total
error(['error - zeta_unsort vector has ',num2str(length(zeta_unsort)),' entries instead of ', ...
num2str(num_modes_total)]);
end
% define dc gains, 31 is head 1, 32 is head 0
omega2 = (2*pi*freqvec)'.^2; % convert to radians and square
% define frequency range for frequency response
freqlo = 501;
freqhi = 25000;
flo=log10(freqlo) ;
fhi=log10(freqhi) ;
f=logspace(flo,fhi,300) ;
frad=f*2*pi ;
% calculate dc gains if uniform damping, peak gains if non-uniform
if zeta_type == 0 % dc gain
gain_h0_coil = abs([force_coil(1)*xn(32,1)/frad(1) ...
force_coil(2:num_modes_total).*xn(32,2:num_modes_total) ...
./omega2(2:num_modes_total)]);
gain_h1_coil = abs([force_coil(1)*xn(31,1)/frad(1) ...
force_coil(2:num_modes_total).*xn(31,2:num_modes_total) ...
./omega2(2:num_modes_total)]);
gain_h0_piezo = abs([force_piezo(1)*xn(31,1)/frad(1) ...
force_piezo(2:num_modes_total).*xn(32,2:num_modes_total) ...
./omega2(2:num_modes_total)]);
gain_h1_piezo = abs([force_piezo(1)*xn(31,1)/frad(1) ...
force_piezo(2:num_modes_total).*xn(31,2:num_modes_total) ...
./omega2(2:num_modes_total)]);
elseif zeta_type == 1 % peak gain
gain_h0_coil = abs([force_coil(1)*xn(32,1)/frad(1) ...
force_coil(2:num_modes_total).*xn(32,2:num_modes_total) ...
./((2*zeta_unsort(2:num_modes_total))'.*omega2(2:num_modes_total))]);
gain_h1_coil = abs([force_coil(1)*xn(31,1)/frad(1) ...
force_coil(2:num_modes_total).*xn(31,2:num_modes_total) ...
./((2*zeta_unsort(2:num_modes_total))'.*omega2(2:num_modes_total))]);
gain_h0_piezo = abs([force_piezo(1)*xn(31,1)/frad(1) ...
force_piezo(2:num_modes_total).*xn(32,2:num_modes_total) ...
./((2*zeta_unsort(2:num_modes_total))'.*omega2(2:num_modes_total))]);
gain_h1_piezo = abs([force_piezo(1)*xn(31,1)/frad(1) ...
force_piezo(2:num_modes_total).*xn(31,2:num_modes_total) ...
./((2*zeta_unsort(2:num_modes_total))'.*omega2(2:num_modes_total))]);
end
% sort gains, keeping track of original and new indices so can rearrange
% eigenvalues and eigenvectors
[gain_h0_coil_sort,index_h0_coil_sort] = sort(gain_h0_coil);
[gain_h1_coil_sort,index_h1_coil_sort] = sort(gain_h1_coil);
[gain_h0_piezo_sort,index_h0_piezo_sort] = sort(gain_h0_piezo);
[gain_h1_piezo_sort,index_h1_piezo_sort] = sort(gain_h1_piezo);
gain_h0_coil_sort = fliplr(gain_h0_coil_sort); % max to min
gain_h1_coil_sort = fliplr(gain_h1_coil_sort); % max to min
gain_h0_piezo_sort = fliplr(gain_h0_piezo_sort); % max to min
gain_h1_piezo_sort = fliplr(gain_h1_piezo_sort); % max to min
index_h0_coil_sort = fliplr(index_h0_coil_sort) % max to min indices
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -