应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 程序语言 » MATLAB/Mathematica » MATLAB编程求解湍流k-e方程时，总是遇见‘下标索引必须为正整数类型或逻辑类型’错误
21/1返回列表
查看: 1727  |  回复: 1
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

且听风吟R

金虫 (小有名气)

[求助] 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=1lthl-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=1lthr-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=1lthb-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=1ltht-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速度错了，不知道是不是，求大神解惑？

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

» 本帖附件资源列表

  • 欢迎监督和反馈：小木虫仅提供交流平台，不对该内容负责。
    本内容由用户自主发布，如果其内容涉及到知识产权问题，其责任在于用户本人，如对版权有异议，请联系邮箱：xiaomuchong@tal.com
  • 附件 1 : lianxi2.m
    • 2017-02-28 20:37:45, 10.82 K
  • 附件 2 : mesh2.m
    • 2017-02-28 20:38:14, 630 bytes
  • 附件 3 : domain.m
    • 2017-02-28 20:38:28, 195 bytes
  • 附件 4 : bubble.m
    • 2017-02-28 20:38:54, 274 bytes
  • 附件 5 : egauss.m
    • 2017-02-28 20:38:54, 1.69 K
  • 附件 6 : gauss2.m
    • 2017-02-28 20:38:55, 1.32 K
  • 附件 7 : result_of.mat
    • 2017-02-28 20:38:55, 6.2 K
  • 附件 8 : shape_1.m
    • 2017-02-28 20:38:55, 1.14 K
  • 附件 9 : weight.m
    • 2017-02-28 20:38:55, 4.12 K

» 猜你喜欢
1楼 2017-02-28 20:34:40
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

且听风吟R

金虫 (小有名气)
下面是几个子函数，这个程序应该是个错误的，现在目的是先把“下标索引必须为正整数类型或逻辑类型”这个错误解决了
赞一下 回复此楼
2楼2017-02-28 20:40:47
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 且听风吟R 的主题更新
21/1返回列表
普通表情 龙 兔 虎 猫
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 初试 317 +6 半拉月丙 2026-03-20 6/300 2026-03-21 18:27 by 学员8dgXkO
[考研] 一志愿重庆大学085700资源与环境总分308求调剂 +7 墨墨漠 2026-03-20 7/350 2026-03-21 16:36 by barlinike
[考研] 22 350 本科985求调剂，求老登收留 +3 李轶男003 2026-03-20 3/150 2026-03-21 13:28 by 搏击518
[考研] 330求调剂0854 +3 assdll 2026-03-21 3/150 2026-03-21 13:01 by 搏击518
[考研] 307求调剂 +3 wyyyqx 2026-03-17 3/150 2026-03-21 03:20 by JourneyLucky
[考研] 初始318分求调剂（有工作经验） +3 1911236844 2026-03-17 3/150 2026-03-21 02:33 by JourneyLucky
[考研] 324分 085600材料化工求调剂 +4 llllkkkhh 2026-03-18 4/200 2026-03-21 01:24 by JourneyLucky
[考研] 一志愿武理材料305分求调剂 +6 想上岸的鲤鱼 2026-03-18 7/350 2026-03-21 01:03 by JourneyLucky
[考研] 274求调剂 +10 S.H1 2026-03-18 10/500 2026-03-20 23:51 by JourneyLucky
[考研] 290求调剂 +7 ^O^乜 2026-03-19 7/350 2026-03-20 21:43 by JourneyLucky
[考研] 260求调剂 +3 朱芷琳 2026-03-20 3/150 2026-03-20 20:35 by 学员8dgXkO
[考研] 319求调剂 +3 小力气珂珂 2026-03-20 3/150 2026-03-20 19:47 by JourneyLucky
[考研] 一志愿吉林大学材料学硕321求调剂 +11 Ymlll 2026-03-18 15/750 2026-03-20 19:40 by 丁丁*
[考博] 申博26年 +3 八6八68 2026-03-19 3/150 2026-03-19 19:43 by nxgogo
[考研] 0854，计算机类招收调剂 +3 胡辣汤放糖 2026-03-15 6/300 2026-03-18 12:09 by 上岸上岸……..
[考研] 0703化学调剂 +3 妮妮ninicgb 2026-03-17 3/150 2026-03-18 10:29 by macy2011
[考研] 283求调剂 +3 听风就是雨； 2026-03-16 3/150 2026-03-17 07:41 by 热情沙漠
[考研] 11408 一志愿西电，277分求调剂 +3 zhouzhen654 2026-03-16 3/150 2026-03-17 07:03 by laoshidan
[考研] 326求调剂 +4 诺贝尔化学奖觊� 2026-03-15 7/350 2026-03-16 17:11 by 诺贝尔化学奖觊�
[考研] 304求调剂 +5 素年祭语 2026-03-15 5/250 2026-03-16 17:00 by 我的船我的海
信息提示
关闭
请填处理意见
关闭