24小时热门版块排行榜

返回列表

且听风吟R

金虫 (小有名气)

应助: 0 (幼儿园)
金币: 636.6
散金: 65
红花: 2
帖子: 205
在线: 84.1小时
虫号: 3198539
注册: 2014-05-11
性别: GG
专业: 流体力学

[求助] MATLAB编程求解湍流k-e方程时，总是遇见‘下标索引必须为正整数类型或逻辑类型’错误

%%%%%%%%%%%%%%%参数设置%%%%%%%%%%%%%
tic
clear
clc
%1.模型参数
LD=1;H=0.04;
mu=0.01;    %动力粘度
rou=1100; %密度
alpha=1e5; %罚函数
alph1=1e11;
sgmk=1.0;
sgms=1.22;
Csgm1=1.44;
Csgm2=1.92;
Cmu=0.09;
%2.边界参数
ul=1;vl=0; %左边界速度
ub=0;vb=0; %下边界
ur=0;vr=0; %右边界
ut=0;vt=0; %上边界

Re=rou*ul*H/mu; %雷诺数
l=0.07*H;       %湍流尺度
I=0.16*Re^(-1/8); %湍流强度
k=1.5*(ul*I)^2; %入口湍动能
s=(Cmu^(3/4))*(k^1.5)/l;  %湍流耗散率

%%%%%%%%%%%%%%输入节点数据%%%%%%%%%%%
ndivlq=20;ndivwq=20;
[x,conn1,conn2,numnod]=mesh3(LD,H,ndivlq,ndivwq);

ndivlq=20;ndivwq=20;       %划背景网格
[xc,conn,numcell,numq]=mesh2(LD,H,ndivlq,ndivwq);

dmax=1.4;                   %放大系数
xspac=LD/ndivlq;          %网格的宽度
yspac=H/ndivwq;             %网格的高度
dm(1,1:numnod)=dmax*xspac*ones(1,numnod);  %x方向节点影响域大小
dm(2,1:numnod)=dmax*yspac*ones(1,numnod);  %y方向节点影响域大小

%%%%%%%%%%%%设置高斯点%%%%%%
%1.湍流场的高斯点
quado=4;             %一维方向上采用的高斯积分点的个数
[gauss]=gauss2(quado);  %函数调用：生成一维方向上的规格化的高斯点
numq2=numcell*quado^2;  %numq2为高斯点总数
gs=zeros(4,numq2);
[gs]=egauss(xc,conn,gauss,numcell);
% plot(gs(1,

,gs(2,

,'.b');
% hold on
% plot(x(1,

,x(2,

,'ro') ;

%%%%%%%%%%%%%%%%%%%%%%%%定义边界节点%%%%%%%%%%%%
ind1=0;ind2=0;ind3=0;ind4=0;
for j=1:numnod
if(x(1,j)==0.0)
      ind1=ind1+1;
      nl(1,ind1)=x(1,j);    %%%左边界
      nl(2,ind1)=x(2,j);
end
if(x(1,j)==LD)
      ind2=ind2+1;
      nr(1,ind2)=x(1,j); %%%右边界
      nr(2,ind2)=x(2,j);
end
if(x(2,j)==0)
      ind3=ind3+1;
      nb(1,ind3)=x(1,j); %%%下边界
      nb(2,ind3)=x(2,j);
end
if(x(2,j)==H)
      ind4=ind4+1;       %%%上边界
      nt(1,ind4)=x(1,j);
      nt(2,ind4)=x(2,j);
end
end

%%===各个边界上节点的个数===%%
lthl=length(nl); %左边界上节点的个数
lthb=length(nb); %下边界上节点的个数
ltht=length(nt); %上边界上节点的个数
lthr=length(nr); %右边界上节点的个数

%%对各边界上的元素进行排序
nl(2,

=bubble(nl(2,

);
nr(2,

=bubble(nr(2,

);
nt(1,

=bubble(nt(1,

);
nb(1,

=bubble(nb(1,

);
%这里还是要进行节点的重新排布，否则高斯点排布就会出错

%设置边界上高斯点的坐标
quado=4;
ind=0;
gauss=gauss2(quado);
for i=1

lthl-1)
ycen=(nl(2,i+1)+nl(2,i))/2;
jcob=abs((nl(2,i+1)-nl(2,i))/2);
for j=1:quado
      mark(j)=ycen-gauss(1,j)*jcob;
      ind=ind+1;
      gsl(1,ind)=nl(1,i);
      gsl(2,ind)=mark(j);
      gsl(3,ind)=gauss(2,j);
      gsl(4,ind)=jcob;
end
end

quado=4;
ind=0;                            %r边界上的高斯点
gauss=gauss2(quado);
for i=1

lthr-1)
ycen=(nr(2,i+1)+nr(2,i))/2;
jcob=abs((nr(2,i+1)-nr(2,i))/2);
for j=1:quado
      mark(j)=ycen-gauss(1,j)*jcob;
      ind=ind+1;
      gsr(1,ind)=nr(1,i);
      gsr(2,ind)=mark(j);
      gsr(3,ind)=gauss(2,j);
      gsr(4,ind)=jcob;
end
end

quado=4;
ind=0;                            %b边界上的高斯点
gauss=gauss2(quado);
for i=1

lthb-1)
ycen=(nb(1,i+1)+nb(1,i))/2;
jcob=abs((nb(1,i+1)-nb(1,i))/2);
for j=1:quado
      mark(j)=ycen-gauss(1,j)*jcob;
      ind=ind+1;
      gsb(1,ind)=mark(j);
      gsb(2,ind)=nb(2,i);
      gsb(3,ind)=gauss(2,j);
      gsb(4,ind)=jcob;
end
end

quado=4;
ind=0;                            %t边界上的高斯点
gauss=gauss2(quado);
for i=1

ltht-1)
ycen=(nt(1,i+1)+nt(1,i))/2;
jcob=abs((nt(1,i+1)-nt(1,i))/2);
for j=1:quado
      mark(j)=ycen-gauss(1,j)*jcob;
      ind=ind+1;
      gst(1,ind)=mark(j);
      gst(2,ind)=nt(2,i);
      gst(3,ind)=gauss(2,j);
      gst(4,ind)=jcob;
end
end

%%%%%%%%%%%%%%%边界积分%%%%%%%%%%%%%%%%%
%1.湍流边界积分
fks=zeros(numnod*2,1);fkspm=zeros(numnod*2,1);
%沿左边界L积分
klk=sparse(numnod,numnod);kls=sparse(numnod,numnod);
ind=0;
for gkl=gsl
ind=ind+1;
gpos=gkl(1:2);
weight=gkl(3);
jac=gkl(4);
v=domain(gpos,x,dm,numnod);
L=length(v);
en=zeros(1,L);vn=zeros(1,L);
flk=zeros(1,L);fls=zeros(1,L);
[phi,dphix,dphiy]=shape_1(gpos,dmax,x,v,dm);
for i=1:L
      en(i)=v(i);vn(i)=v(i);
      flk(i)=k*phi(i);
      fls(i)=s*phi(i);
end
klk(en,en)=klk(en,en)+sparse((weight*jac)*(phi'*phi));
kls(en,en)=kls(en,en)+sparse((weight*jac)*(phi'*phi));
fkspm(en)=fkspm(en)+jac*weight*flk';
fkspm(en+numnod)=fkspm(en+numnod)+jac*weight*fls';
end
%沿下边界B积分
for gksb=gsb
gpos=gksb(1:2);
weight=gksb(3);
jac=gksb(4);
v=domain(gpos,x,dm,numnod);
L=length(v);
en=zeros(1,L);
fkb=zeros(1,L);fsb=zeros(1,L);
Kb=0;Sb=0;
[phi,dphix,dphiy]=shape_1(gpos,dmax,x,v,dm);
for i=1:L
      en(i)=v(i);
      fkb(i)=Kb*phi(i);
      fsb(i)=Sb*phi(i);
end
fks(en)=fks(en)+jac*weight*fkb';
fks(en+numnod)=fks(en+numnod)+jac*weight*fsb';
end
%沿上边界T积分
for gkst=gst
gpos=gkst(1:2);
weight=gkst(3);
jac=gkst(4);
v=domain(gpos,x,dm,numnod);
L=length(v);
en=zeros(1,L);
fkt=zeros(1,L);fst=zeros(1,L);
Kt=0;St=0;
[phi,dphix,dphiy]=shape_1(gpos,dmax,x,v,dm);
for i=1:L
      en(i)=v(i);
      fkb(i)=Kt*phi(i);
      fsb(i)=St*phi(i);
end
fks(en)=fks(en)+jac*weight*fkt';
fks(en+numnod)=fks(en+numnod)+jac*weight*fst';
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 非牛顿迭代 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1、计算初始黏度(viscosity)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 迭代 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1.进行初始赋值
load ('result_of.mat'); %调用相同边界条件下NS方程解速度分布结果ux，uy
kkk=ones(1,numnod);    %湍动能初始值
kks=ones(1,numnod);    %耗散率初始值
vis=ones(1,numnod);

vis_k_1=vis;

kkk_k_1=kkk;
kks_k_1=kks;

%%%%%%%%%%%%%% 主程序迭代开始 %%%%%%%%%%%%%%%%
%1.迭代初始条件
norm_vis=1;
norm_k=1;
norm_s=1;
times=0;

fprintf('  现在开始计算，请耐心等待  \n')
while ((norm_vis>1e-3||norm_vis>1e-3||norm_k>1e-3||norm_s>1e-3)&&times<1000)
%迭代赋值:将k+1时的结果赋值给k时的值
kkk_k=kkk_k_1;
kks_k=kks_k_1;
vis_k=vis_k_1;

%%%%%%%%%%%%%%%%%%%%% 对高斯点循环，组装K矩阵 %%%%%%%%%%%%%%%%%%%%%%%%%%
% 1、更新对流项子块
% ！！开始前需计算温度场上节点速度,本程序速度场与温度场节点重合,所以不需重新计算，直接利用速度场上节点的速度！
%%%%%%%%%%%%%%%%%%%%%%%%%%%计算k-e方程%%%%%%%%%%%%%%%%%%
%1.计算湍流方程扩散项K矩阵
kcl=sparse(numnod,numnod);
kcs=sparse(numnod,numnod);
for ggk=gs
gpos=ggk(1:2);
weight=ggk(3);
jac=ggk(4);
vk=domain(gpos,x,dm,numnod);
Lk=length(vk);
en=zeros(1,Lk);vn=zeros(1,Lk);
[phi,dphix,dphiy]=shape_1(gpos,dmax,x,vk,dm);
      for i=1:Lk
         en(i)=vk(i);
      end
      kcl(en,en)=kcl(en,en)+sparse((vis_k(1,i)/sgmk)*(weight*jac)*(dphix'*dphix+dphiy'*dphiy));
      kcs(en,en)=kcs(en,en)+sparse((vis_k(1,i)/sgmk)*(weight*jac)*(dphix'*dphix+dphiy'*dphiy));
end

%2.计算湍流方程对流项K矩阵,
%开始前需计算湍流场上节点速度,本程序速度场与湍流场节点重合,所以不需重新计算，直接利用速度场上节点的速度！
kdl=sparse(numnod,numnod);
kds=sparse(numnod,numnod);
for ggs=gs
      gpos=ggs(1:2);
      weight=ggs(3);
      jac=ggs(4);
      v=domain(gpos,x,dm,numnod);
      L=length(v);
      en=zeros(1,L);
      [phi,dphix,dphiy]=shape_1(gpos,dmax,x,v,dm);
      for i=L
         en(i)=v(i);
      end
      kdl(en,en)=kdl(en,en)+sparse((weight*jac)*(phi*phi'*ux_k_1*dphix'+phi*phi'*uy_k_1*dphiy'));
      kds(en,en)=kds(en,en)+sparse((weight*jac)*(phi*phi'*ux_k_1*dphix'+phi*phi'*uy_k_1*dphiy'));
end

%3.计算耗散项K矩阵
%开始前需计算湍流场上节点速度,本程序速度场与湍流场节点重合,所以不需重新计算，直接利用速度场上节点的速度！
krk=sqarse(numnod,numnod);
krs=sqarse(numnod,numnod);
for ggu=gs
      gpos=ggu(1:2);
      weight=ggu(3);
      jac=ggu(4);
      v=domain(gpos,x,dm,numnod);
      L=length(v);
      en=zeros(1,L);
      [phi,dphix,dphiy]=shape_1(gpos,dmax,x,v,dm);
      for i=L
         en(i)=v(i);
      end
      krk(en,en)=krk(en,en)+sparse((weight*jac)*vis_k(1,i)*phi*(2*(dphix'*ux_k_1*dphix'*ux_k_1)+2*(dphiy'*uy_k_1*dphiy'*uy_k_1)+(dphiy'*ux_k_1*dphiy'*ux_k_1+2*dphiy'*ux_k_1*dphix'*uy_k_1+dphix'*uy_k_1*dphix'*uy_k_1)));
      krs(en,en)=krs(en,en)+sparse((weight*jac)*Csgm1*yinzi_k(1,i)*vis_k(1,i)*phi*(2*(dphix'*ux_k_1*dphix'*ux_k_1)+2*(dphiy'*uy_k_1*dphiy'*uy_k_1)+(dphiy'*ux_k_1*dphiy'*ux_k_1+2*dphiy'*ux_k_1*dphix'*uy_k_1+dphix'*uy_k_1*dphix'*uy_k_1)));
end

%4.计算湍流M矩阵
Mks=sparse(numnod,numnod);
Msk=sparse(numnod,numnod);
for ggu=gs
      gpos=ggu(1:2);
      weight=ggu(3);
      jac=ggu(4);
      v=domain(gpos,x,dm,numnod);
      L=length(v);
      en=zeros(1,L);
      [phi,dphix,dphiy]=shape_1(gpos,dmax,x,v,dm);
      for i=L
         en(i)=v(i);
      end
      Mks(en,en)=Mks(en,en)+sparse((weight*jac)*(phi'*phi));
      Msk(en,en)=Msk(en,en)+sparse((weight*jac)*Csgm2*yinzi_k(1,i)*(phi'*phi));
end
%5.构建k-e总体方程
k_ks=[kdl-kcl,Mks;zeros(numnod,numnod),kds-kcs+Msk];
k_s=[krk;krs];
k_b=[klk;kls];
f_ks=fks+k_s+k_b+aplh*fkspm;
s=k_ks\f_ks;
ks=d(1:numnod*2);
for i=1:numnod
      K(i)=ks(i);
      E(i)=ks(i+numnod);
end

%6.更新粘度
vis_k_1=zeros(1,numnod);          % 初始化vis
yinzi_k_1=zeros(1,numnod);
for gg=gs
      gpos=gg(1:2);
      weight=gg(3);
      jac=gg(4);
      v=domain(gpos,x,dm,numnod);
      L=length(v);
      en=zeros(1,L);
      kk=zeros(1,L);
      ks=zeros(1,L);
      [phi,dphix,dphiy]=shape_1(gpos,dmax,x,v,dm);
      for i=1:L
         en(i)=v(i);
         ke=phi(i)*K(i);
         kk=(phi(i)'*phi(i))*(K(i))^2;
         ks=(phi(i))*E(i);
      end
      vis_k_1(en)=vis_k_1(en)+(kk/ks);
      yinzi_k_1(en)=yinzi_k_1(en)+(ks/ke);
end
   if norm(kkk_k_1-kkk_k)<1e-3
      norm_kkk=0;
else
      norm_kkk=norm(kkk_k_1-kkk_k)/norm(kkk_k);
end
if norm(kks_k_1-kks_k)<1e-3
      norm_kks=0;
else
      norm_kks=norm(kks_k_1-kks_k)/norm(kks_k);
end
if norm(vis_k_1-vis_k)<1e-3
      norm_vis=0;
else
      norm_vis=norm(vis_k_1-vis_k)/norm(vis_k);
end
if norm(yinzi_k_1-yinzi_k)<1e-3
      norm_yinzi=0;
else
      norm_yinzi=norm(yinzi_k_1-yinzi_k)/norm(yinzi_k);
end

clear k ku11 ku12 ku21 ku22 kp11 kp12 kp21 kp22
% 7.输出迭代结果
times=times+1;
fprintf('time= %d  && norm_kkk= %6.9f && norm_kks= %6.9f && norm_vis= %6.9f && norm_yinzi= %6.9f \n',times,norm_kkk,norm_kks,norm_vis,norm_yinzi)
end

前面的NS方程检查过没有问题，但是在解k-e方程的时候，在kdl矩阵那出现了“下标索引必须为正整数类型或逻辑类型”错误，我觉得是调用ux-k-1速度错了，不知道是不是，求大神解惑？

后面应该还有很多问题，先把这个给解决再说

回复此楼