%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Proper Orthogonal Decomposition
%参考文献:JFM 2007 Vol 583,pp 199-227
%A Turbulent Jet in Crossflow Analysed with Proper Orthogonal Decomposition
%K. E. Meyer et al.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
fprintf('=============================================================\n');
idim=901;
jdim=1;
kdim=961;
nstart=605;
nend=624;
%nend=1746;
N=nend-nstart+1;
M=idim*jdim*kdim;
dt=1.0/N;
nvar=2;
u=zeros(M,N);
% dlmread('filename',delimiter,row,col) 将以ASCII码分隔的数值数据读入到矩阵中
% 其中参数delimiter用于指定文件中的分隔符,起始行为row,起始列为col
fid=fopen('.\2Doutyz\grid.bin','rb');
[grid,num]=fread(fid,inf,'float');
fprintf('读取时间序列数据:\n');
fclose(fid);
ncount=0;
for n=nstart:nend
ncount=ncount+1
ch=num2str(n,'%4.4d');
% 带变化字符的fid需要加入方括号[]。
fid=fopen(['.\2Doutyz\2ddata' ch '.bin'],'rb');
[Q,num]=fread(fid,inf,'float');
for k=1:kdim
for i=1:idim
num1=i+(k-1)*idim+1;
num2=num1+M*(nvar-1);
u(num1-1,ncount)=Q(num2);
end
end
fclose(fid);
clear Q;
end
uave=zeros(M,1);
for n=1:ncount
uave(1:M,1)=uave(1:M,1)+u(1:M,n);
end
uave(1:M,1)=uave(1:M,1)/ncount;
for n=1:ncount
u(1:M,n)=-uave(1:M,1)+u(1:M,n);
end
clear uave;
fid=fopen('.\test.dat','wt');
fprintf(fid,'variables = "x", "r", "u"\n');
fprintf(fid,'ZONE I= %6.5d , K= %6.5d F=POINT\n',idim,kdim);
for k=1:M
fprintf(fid, '%10.7f %10.7f %10.7f \n',grid(k+1),grid(k+1+M),u(k,ncount));
end
fclose(fid);
A=u(1:M,1:N)'*u(1:M,1:N);
%end
fprintf('求特征值和特征向量:\n');
%V是特征矩阵,其每一列是对应的特征向量,D是升序排列的特征值
[V,D]=eig(A);
fprintf('按特征值降序排列:\n');
%调整特征向量V和特征值D,对于对称矩阵求的特征值是自动排序(由小到大),我们需要
%按降序排列,
Vtemp=V;
Dtemp=D;
for n=1:N
V(1:N,n)=Vtemp(1:N,N-n+1);
D(n,n)=Dtemp(N-n+1,N-n+1);
end
clear Vtemp;
clear Dtemp;
% 验证
tem1=zeros(N,1);
tem2=zeros(N,1);
n=1;
tem1(1:N,1)=V(1:N,n)*D(n,n)
tem2(1:N,1)=A(1:N,1:N)*V(1:N,n)
clear A;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%输出特征值
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('输出特征值:\n');
d=diag(D);
d=d/sum(d); %特征值的归一化处理!!
fid=fopen('.\PODresult\POD_eigne.dat','w');
fprintf(fid,'variables = Mode, lamda\n');
for n=1:N
fprintf(fid, '%12.8d %12.8f\n',n,d(n));
end
fclose(fid);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('求模态:\n');
Phi=zeros(M,N);
for n=1:N
Phi(1:M,n)=u(1:M,1:N)*V(1:N,n); %特征向量上的投影
Phi(1:M,n)=Phi(1:M,n)/norm(Phi(1:M,n)); %Phi的每一列为对应的模态
end
fprintf('求系数:\n');
a(N,N)=0.0;
for n=1:N
a(1:N,n)=Phi(1:M,1:N)'*u(1:M,n); %%%得到系数矩阵~
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%输出模态和模态系数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('输出模态及系数:\n');
for n=1:N
fprintf('mode:%2.2d\n',n);
ch = num2str(n,'%2.2d');
%%%%%%%%%%%%%%%%% 模态向量!!
fid=fopen(['.\POD_mode\POD_mode' ch '.plt'],'wt');
fprintf(fid,'variables = "x", "z", "mode"\n');
fprintf(fid,'ZONE I= %6.5d , K= %6.5d F=POINT\n',idim,kdim);
for k=1:M
fprintf(fid, '%10.7f %10.7f %10.7f \n',grid(k+1),grid(k+1+M),Phi(k,n));
end
fclose(fid);
%%%%%%%%%%%%%%%%% 模态系数!!!
fid=fopen(['.\POD_coef\POD_coef' ch '.dat'],'w');
time=1:N;%time为列向量表征1-N模态
timer=time*0.02;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%非常重要!!!!!!!!!!!!!!!!
aa=a';
%%%%%注意,这里讲每一个模态对应的系数,
%%%%%即此JFM文中,a(i,n)*Phi(:,i),要给出与i对应的a,因此需要转置。
%%%%%%%%%%%%%非常重要!!!!!!!!!!!!!!!!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
coeffi=aa(:,n);
coeffi=[time;coeffi'];
fprintf(fid,'%7.5f,%12.8f\n',coeffi);
fclose(fid);
end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%重构流畅
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%
% %%
%
% %%
% Un(M,N)=0;
% Um=zeros(M,3);
% for n=1:N
% fprintf('mode:%2.2d\n',n);
% %Un(1:M,n)=Umean(1:M);
% for i=1:4 %这里可以减少重构的模态
% Un(1:M,n)=Un(1:M,n)+Phi(1:M,i)*a(i,n);
% end
% Um(1:M,1)=Um(1:M,1)+Un(1:M,n);
% Um(1:M,2)=Um(1:M,2)+u(1:M,n)-Un(1:M,n);
% Um(1:M,3)=Um(1:M,3)+u(1:M,n);
% end
%
% Um(1:M,1)=Um(1:M,1)/N;
% Um(1:M,2)=Um(1:M,2)/N;
% Um(1:M,3)=Um(1:M,3)/N;
%
% fprintf('输出重构流场:\n');
% fid=fopen(['.\PODrecon\POD_recon4.dat'],'wt');
% fprintf(fid,'variables = "x", "z" , "u_r", "u_s", "u"\n');
% fprintf(fid,'ZONE I= %6.5d , K= %6.5d F=POINT\n',idim,kdim);
% for k=1:M
% fprintf(fid, '%10.7f %10.7f %10.7f \n',grid(k+1),grid(k+1+M),Um(k,1),Um(k,2),Um(k,3));
% end
% fclose(fid);
%
% %%
%fprintf('验证重构是否正确:\n');
%fid=fopen(['POD_Validation.dat'],'wt');
%for n=1:N
% fprintf('mode:%2.2d\n',n);
% for m=1:M
% fprintf(fid, '%10.7f\n',Un(m,n)-u(m,n));
% end
%end
%fclose(fid);
% 计算模态能量
fprintf('计算模态能量\n');
enn=0;
En=zeros(N);
for n=1:N
for m=1:N
for i=1:M
En(n)=En(n)+(a(n,m)*Phi(i,m))^2;
end
end
enn=enn+En(n);
end
%ch = num2str(n,'%2.2d');
fid=fopen(['.\PODE\POD_En.dat'],'w');
for n=1:N
En(n)=En(n)/enn;
fprintf(fid, '%7.5d,%10.7f\n',n,En(n));
end
fclose(fid)
%
clear
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
没有合适的资源?快使用搜索试试~ 我知道了~
POD.rar_POD_POD analysis_combinationzwe_pod 波氏分析_rodscj
共1个文件
m:1个
1.该资源内容由用户上传,如若侵权请联系客服进行举报
2.虚拟产品一经售出概不退款(资源遇到问题,请及时私信上传者)
2.虚拟产品一经售出概不退款(资源遇到问题,请及时私信上传者)
版权申诉
0 下载量 123 浏览量
2022-07-14
05:06:03
上传
评论
收藏 2KB RAR 举报
温馨提示
可以用于POD分析,得到不同模态以及重建流场
资源推荐
资源详情
资源评论
收起资源包目录
POD.rar (1个子文件)
POD.m 6KB
共 1 条
- 1
资源评论
邓凌佳
- 粉丝: 65
- 资源: 1万+
上传资源 快速赚钱
- 我的内容管理 展开
- 我的资源 快来上传第一个资源
- 我的收益 登录查看自己的收益
- 我的积分 登录查看自己的积分
- 我的C币 登录后查看C币余额
- 我的收藏
- 我的下载
- 下载帮助
安全验证
文档复制为VIP权益,开通VIP直接复制
信息提交成功