%%%%%%%%%%%%%%%参数设置%%%%%%%%%%%%%
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)&×<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速度错了,不知道是不是,求大神解惑?
后面应该还有很多问题,先把这个给解决再说