% PCKFCURVES Traces the dispersion curves of a waveguide in the K-F space
%
% A lot of work has been put into optimizing this routine because tracing
% the curves takes a lot of time, but remember that "you're bound to be
% unhappy if you try to optimize everything" (Donald Knuth)
%
% PCDISP: Pochhammer-Chree dispersion in cylindrical waveguides
% Written by Fernando Seco Granja (fernando.seco@csic.es)
% Consejo Superior de Investigaciones Cient�ficas (CSIC)
% http://www.car.upm-csic.es/lopsi/people/fernando.seco/pcdisp/
%
% Last update: 3-10-2011
clear,clc
% ********** PARTS OF THE SPECTRUM TO BE PLOTTED ********
flag_realroots=1;
flag_imagroots=1;
flag_complexroots=1;
% ********** WAVEGUIDE PARAMETERS **********
pcwaveguide
% Circumferential order to be plotted (n = -1 for torsional modes T(0,m),
% n = 0 for longitudinal L(0,m), and n >= 1 for flexural, F(n,m))
n = 0;
% *********** SOLVER PARAMETERS *********
% Solutions will be computed in the 3D space with 0<Re{k},Im{k}<=kmax,
% and 0<f<=fmax
fmax=1000e3; % Maximum frequency (Hz)
kmax=2000; % Maximum wavenumber (1/m)
% Tolerance error in the determination of the roots, in normalized units
drtol=1e-6;
% Maximum distance between consecutive points along a branch
drmax=.05;
drmin=5*drtol;
dr_factor=1.33; % should be strictly bigger than 1
% Other solver parameters
dthmax=pi/4;
Nroots=20;
% Number of points between consecutive roots tested for possible sign
% changes in the bisection root finder. It should be increased if roots are
% missing because they're very close
npcheck=150;
% ********* HERE WE GO! **********
% Normalized coordinates for 2D plot (x -> k, y -> f)
xmax=kmax/sc_k;
ymax=fmax/sc_freq;
clc
disp('PCKFCURVES running')
disp(' ')
figure(1),clf,hold on
axis(1.01*[-kmax kmax 0 fmax/1e3])
xlabel('-Im\{k\} (1/m) Re\{k\}')
ylabel('f (kHz)')
handle_branch=[];
indbranch=0;
xbranch=[];
ybranch=[];
npbranch=[];
flagbranch=[];
% ****** TRACING OF BRANCHES WITH REAL AND IMAGINARY WAVENUMBERS ******
tvect=[];
if(flag_realroots),
tvect=[1];
end
if(flag_imagroots),
tvect=[tvect j];
end
for t=tvect,
if(t==1),
color='b';xsign=1;name='real';
else,
color='r';xsign=-1;name='imaginary';
end
% First we collect the roots of the frequency equation along the edges
% of the wavenumber and frequency boxes. These are the starting and
% ending points for the branches, and are called "docking" points here
disp(sprintf('Determining docking points for branches with %s wavenumber...',name))
dr=drmax/3;
np=3; % number of pre-computed points at the beginning/end of each branch
numleft=[];numright=[];numtop=[];numbottom=[];numcircle=[];
for ind=1:np,
% Left edge of the k-f space (k=0)
dummy=pcsolvebisection('k-f',t*dr*ind,[4/3*drmax ymax],drtol,-1,npcheck,n,waveguidedescription);
numleft(ind)=length(dummy);
yleft(ind,1:numleft(ind))=dummy;
xleft(ind,1:numleft(ind))=dr*ind;
% Right edge of the k-f space (k=kmax)
dummy=pcsolvebisection('k-f',t*(xmax-dr*(ind-1)),[dr ymax],drtol,-1,npcheck,n,waveguidedescription);
numright(ind)=length(dummy);
yright(ind,1:numright(ind))=dummy;
xright(ind,1:numright(ind))=xmax-dr*(ind-1);
% Top edge of the k-f space (f=fmax)
dummy=pcsolvebisection('f-k',ymax-dr*(ind-1),t*[0 xmax],drtol,-1,npcheck,n,waveguidedescription);
numtop(ind)=length(dummy);
xtop(ind,1:numtop(ind))=dummy;
ytop(ind,1:numtop(ind))=ymax-dr*(ind-1);
% Bottom edge of the k-f space (f=0)
dummy=pcsolvebisection('f-k',dr*ind,t*[4/3*drmax xmax],drtol,-1,npcheck,n,waveguidedescription);
numbottom(ind)=length(dummy);
xbottom(ind,1:numbottom(ind))=dummy;
ybottom(ind,1:numbottom(ind))=dr*ind;
% Quarter circle around the origin (k=0,f=0)
dummy=pcsolvebisection('0-0',t*dr*ind,[0 pi/2],drtol/(ind*dr),-1,npcheck,n,waveguidedescription);
numcircle(ind)=length(dummy);
xcircle(ind,1:numcircle(ind))=dr*ind*cos(dummy);
ycircle(ind,1:numcircle(ind))=dr*ind*sin(dummy);
end
flag_plot_dockpoints=0;
if(flag_plot_dockpoints),
plot(xleft*sc_k*xsign,yleft*sc_freq/1e3,'Marker','o','Color',color)
plot(xright*sc_k*xsign,yright*sc_freq/1e3,'Marker','x','Color',color)
plot(xtop*sc_k*xsign,ytop*sc_freq/1e3,'Marker','+','Color',color)
plot(xbottom*sc_k*xsign,ybottom*sc_freq/1e3,'Marker','*','Color',color)
plot(xcircle*sc_k*xsign,ycircle*sc_freq/1e3,'Marker','s','Color',color)
end
% Checks that the starting points have been correctly computed
if(any(numleft(2:np)~=numleft(1)) | any(numright(2:np)~=numright(1)) | ...
any(numtop(2:np)~=numtop(1)) | any(numbottom(2:np)~=numbottom(1)) | ...
any(numcircle(2:np)~=numcircle(1))),
numleft,numright,numtop,numbottom,numcircle
disp('Problems determining the starting points')
stop
else,
numleft=numleft(1);numright=numright(1);numtop=numtop(1);
numbottom=numbottom(1);numcircle=numcircle(1);
end
% Ocassionally the branch starting at the origin might have been
% computed twice - we have to remove one instance
if(numcircle==1),
slope1=sum(ycircle(1:np,1))/sum(xcircle(1:np,1));
if(slope1<1/3),
if(numbottom>0),
xbottom=xbottom(1:np,2:numbottom);
ybottom=ybottom(1:np,2:numbottom);
numbottom=numbottom-1;
end
end
if(slope1>3),
if(numleft>0),
xleft=xleft(1:np,2:numleft);
yleft=yleft(1:np,2:numleft);
numleft=numleft-1;
end
end
end
% Simple check
if(t==1),
if(numbottom>0 | numleft+numcircle~=numright+numtop),
disp('There was a problem computing the docking points')
numleft,numbottom,numright,numtop,numcircle
stop
end
end
% Arrays with the docking points
xdock=[xleft(1:np,1:numleft) xright(1:np,1:numright) xtop(1:np,1:numtop) xbottom(1:np,1:numbottom) xcircle(1:np,1:numcircle)]';
ydock=[yleft(1:np,1:numleft) yright(1:np,1:numright) ytop(1:np,1:numtop) ybottom(1:np,1:numbottom) ycircle(1:np,1:numcircle)]';
Ndock=size(xdock,1);
% Sort them
[dummy,indord]=sort(xdock(:,1).^2+ydock(:,1).^2);
xdock=xdock(indord,:);ydock=ydock(indord,:);
% We plot them and get a handle for them in the figure
figure(1),
handle_dock=plot(xdock(:,1)*xsign*sc_k,ydock(:,1)*sc_freq/1e3,...
'Marker','o','Color',color,'LineStyle','none');
% Total number of branches
Nbranches=Ndock/2;
if(floor(Nbranches)~=Nbranches),
disp('The number of docking points should be even!')
stop
else,
disp(sprintf('Tracing %g branches with %s wavenumber',Nbranches,name))
end
% Loop over all branches
for indbranchtmp=1:Nbranches,
numeval=0;
xp=xdock(1,1:3);
yp=ydock(1,1:3);
% We eliminate the starting points from the arrays
xdock=xdock(2:end,:);ydock=ydock(2:end,:);
Ndock=Ndock-1;
indbranch=indbranch+1;
flagbranch(indbranch)=t;
figure(1),handle_branch(indbranch)=plot(xp*xsign*sc_k,yp*sc_freq/1e3,'Marker','.','Color',color);
flag_branch_complete=0;
indp=3;
drlast=[norm([xp(3)-xp(2) yp(3)-yp(2)])];
drlista=[drlast*dr_factor drlast drlast/dr_factor];
% When tracing the branch we keep track of the sign of the
% determinant at the left side of the curve
没有合适的资源?快使用搜索试试~ 我知道了~
温馨提示
关于波导杆、管等器件的频散特性曲线仿真,可以根据窗口,设置对应物理属性,密度,直径、孔直径、泊松比等,进行相关的频散曲线仿真,可以有效代替一些仿真软件。 参数设置 material-1,核心层材料,默认钢 弹性模量 泊松比 密度 material 2,包裹层材料,默认水泥砂浆 材料编号 1层材料1 2层结构1 2 半径向量 第一个为空心半径,若实心则为0 其余为各层外半径 如12 mm钢筋,输入0 6 频率范围:[min step max],以kHz为单位 模态选择: -1:T模态 0:L模态 1:F模态 若需要求解多种,则输入: -1 0 1。 半径参数约定 single layer radius-5mm solid bar, radiilayer = [0 5] embeded solid bar, radiilayer = [0 5 8] infinitely embeded solid bar, radiilayer = [0 5 inf] hollow tube ,radiilayer = [5 10] 给定必要参数后,程序执行计算,可以输出频散曲线。
资源详情
资源评论
资源推荐
收起资源包目录
1612_PCdisp_Update.rar (53个子文件)
1612_PCdisp_Update
pcplotmatdet2Dvar.m 2KB
pcsignalpropagation.m 5KB
pcwaveform.m 11KB
pcplotexcitation.m 5KB
signalbank.m 4KB
pcdisp2.m 9KB
pcwaveguide.m 5KB
pcplotmatdet2D.m 2KB
pcspeeds.m 2KB
pckfcurves.m 33KB
tools.m 15KB
PCdisp_源程序
pcplotmatdet2Dvar.m 2KB
pcsignalpropagation.m 5KB
pcwaveform.m 11KB
pcplotexcitation.m 5KB
signalbank.m 4KB
pcwaveguide.m 5KB
pcplotmatdet2D.m 2KB
pcspeeds.m 2KB
pckfcurves.m 33KB
pcdisp.m 9KB
pcmodalanalysis.m 13KB
intsimpson.m 764B
pcextsurfacestress.m 1KB
pcextvolumforce.m 798B
pcplotmatdet1D.m 4KB
intsimpson2.m 841B
pcviewmatdet.m 4KB
pcmat.m 16KB
pcmatdet.m 20KB
pcsolverandom.m 5KB
pcmodename.m 1KB
pcsolvebisection.m 5KB
pccreatefrequencyvect.m 607B
pcorthogonalcheck.m 11KB
pcdisp.m 9KB
pcmodalanalysis.m 13KB
intsimpson.m 764B
pcwaveform2.m 11KB
pcextsurfacestress.m 1KB
pcwaveguide2.m 6KB
pcextvolumforce.m 798B
pcplotmatdet1D.m 4KB
intsimpson2.m 841B
pcviewmatdet.m 4KB
pcmat.m 16KB
pcmatdet.m 20KB
pcsolverandom.m 5KB
pcmodename.m 1KB
pcsolvebisection.m 5KB
pccreatefrequencyvect.m 607B
pcwaveform3.m 8KB
pcorthogonalcheck.m 11KB
共 53 条
- 1
weixin_52874170
- 粉丝: 0
- 资源: 1
上传资源 快速赚钱
- 我的内容管理 展开
- 我的资源 快来上传第一个资源
- 我的收益 登录查看自己的收益
- 我的积分 登录查看自己的积分
- 我的C币 登录后查看C币余额
- 我的收藏
- 我的下载
- 下载帮助
安全验证
文档复制为VIP权益,开通VIP直接复制
信息提交成功
评论0