%%%%%%%%%%%%%%%参数设置%%%%%%%%%%%%%

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速度错了，不知道是不是，求大神解惑？

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