24小时热门版块排行榜

返回列表

当前主题已经存档。

【悬赏金币】回答本帖问题，作者hotoil将赠送您 60 个金币

hotoil

铜虫 (正式写手)

应助: 1 (幼儿园)
金币: 894.5
散金: 153
红花: 1
帖子: 404
在线: 199.6小时
虫号: 412432
注册: 2007-06-25
性别: GG
专业: 多相流热物理学

[求助] 【求助：matlab程序提速】+【有效期至2009年01月03日】，

本人编写了采用TVD格式的Lax-Friedrichs算法求解偏微分方程的matlab程序，请高人帮助提高该程序的运行速度！有问题联系qq350204939
syms q1 q2 q3;
denL=1000;
D=0.146;
A=pi/4*D^2;
qIG=0;
qIL=0;
ang=0;
z=1;
R=100000/300;%J/(kg*K),空气;R=Rm*1000/M_G=8.3143*1000/M_G.
T=273.15+5;
XL=10000;
tstep=0.01;
xstep=200;
m=XL/xstep;
t_N=2400;
fw=0.02;
g=9.8;
dpdx=2/D*fw*505*(0.04*4/pi/D^2)^2;
Q=[q1;q2;q3];%其中 y=u*h;
F=[q1*q3/(q1+q2);q2*q3/(q1+q2);q3^2/(q1+q2)+z*R*T*q1/(1-q2/denL)];
S=[q1/(1-q2/denL)*qIG/A;denL*qIL/A; -2/D*fw*q3^2/(q1+q2)-(q1+q2)*g*sin(ang)];
dFdQ=jacobian(F,Q);
lamda=eig(dFdQ);
lamdaJ=ones(3,m);
FL=ones(3,1);
FLM=ones(3,1);
FR=ones(3,1);
FRM=ones(3,1);
F_TVDLFJ=ones(3,1);
F_TVDLFJ1=ones(3,1);
OLR=ones(3,1);
QJnmid=ones(3,1);
QJ1nmid=ones(3,1);
QJn1=ones(3,1);
QL=ones(3,1);
QQJ=ones(3,1);
QQJ1=ones(3,1);
QR=ones(3,1);
Qtem=ones(3,1);
SJ=ones(3,m+1);
Xtem=ones(3,1);
Ytem=ones(3,1);
vdenGn=ones(1,m+2);
vEln=ones(1,m+2);
vUmn=ones(1,m+2);
vdenGn1=ones(1,m+2);
vEln1=ones(1,m+2);
vUmn1=ones(1,m+2);
Ptem=1;
Ptem1=1;
Ptem2=1;
Ptem3=1;
qq1=1;
qq2=1;
qq3=1;
lamd_max=1;
lamdn_max=1;
n=1;
j=1;
for (j=2:m+1)
Ptem=10*100000+(XL-(j-1.5)*xstep)*dpdx;
vdenGn(j)=Ptem/(z*R*T);
vEln(j)=0.5;
vUmn(j)=0.04*4/pi/D^2;
end
Ptem=10*100000+XL*dpdx;
vdenGn(1)=Ptem/(z*R*T);
vEln(1)=0.5;
vUmn(1)=0.04*4/pi/D^2;
QJ(:,1)=subs(Q,{q1,q2,q3},{vdenGn(1)*(1-vEln(1)),denL*vEln(1),(vdenGn(1)*(1-vEln(1))+denL*vEln(1))*vUmn(1)});
Ptem=10*100000;
vdenGn(m+2)=Ptem/(z*R*T);
vEln(m+2)=0.5;
vUmn(m+2)=0.04*4/pi/D^2;
QJ(:,m+2)=subs(Q,{q1,q2,q3},{vdenGn(m+2)*(1-vEln(m+2)),denL*vEln(m+2),(vdenGn(m+2)*(1-vEln(m+2))+denL*vEln(m+2))*vUmn(m+2)});
vtn=0;
while (vtn lamdn_max=1;
if (vtn+tstep<=1000)
      vEln1(1)=0.5;
      vUmn1(1)=0.04/A;
else
      if(vtn+tstep>1000&vtn+tstep<1010)
         vEln1(1)=1/(2+(n+1-1000/tstep)*tstep/10);
         vUmn1(1)=0.02/A*(2+(n+1-1000/tstep)*tstep/10);;
      else
         vEln1(1)=1/3;
         vUmn1(1)=0.06/A;
      end
end%已知边界条件
vdenGn1(1)=vdenGn(1);
F_TVDLFJ=subs(F,{q1,q2,q3},{vdenGn1(1)*(1-vEln1(1)),denL*vEln1(1),(vdenGn1(1)*(1-vEln1(1))+denL*vEln1(1))*vUmn1(1)});
for (j=2:m)
      qq1=vdenGn(j)*(1-vEln(j));
      qq2=denL*vEln(j);
      qq3=(vdenGn(j)*(1-vEln(j))+denL*vEln(j))*vUmn(j);
      lamdaJ(:,j)=subs(lamda,{q1,q2,q3},{qq1,qq2,qq3});
      QJ(:,j)=subs(Q,{q1,q2,q3},{qq1,qq2,qq3});
      SJ(:,j)=subs(S,{q1,q2,q3},{qq1,qq2,qq3});
      qq1=vdenGn(j+1)*(1-vEln(j+1));
      qq2=denL*vEln(j+1);
      qq3=(vdenGn(j+1)*(1-vEln(j+1))+denL*vEln(j+1))*vUmn(j+1);
      QJ(:,j+1)=subs(Q,{q1,q2,q3},{qq1,qq2,qq3});
      SJ(:,j+1)=subs(S,{q1,q2,q3},{qq1,qq2,qq3});
      if(j>2)%计算QQ(:,j)
         Xtem=QJ(:,j)-QJ(:,j-1);
         Ytem=QJ(:,j+1)-QJ(:,j);
      else%if(j==2)
         Ytem=QJ(:,j+1)-QJ(:,j);
         Xtem=Ytem;%(Xtem+Ytem)/2;
      end
      for(k=1:3)
      QQJ(k)=1/2*(sign(Xtem(k))+sign(Ytem(k)))*min(abs(Xtem(k)),abs(Ytem(k)));
      end
      if(j          Xtem=QJ(:,j+1)-QJ(:,j);
         Ytem=QJ(:,j+2)-QJ(:,j+1);
      else%if(j==m)
         Xtem=QJ(:,j+1)-QJ(:,j);
         Ytem=Xtem;%(Xtem+Ytem)/2;
      end
      for(k=1:3)
      QQJ1(k)=1/2*(sign(Xtem(k))+sign(Ytem(k)))*min(abs(Xtem(k)),abs(Ytem(k)));
      end
      Qtem=QJ(:,j)+1/2*QQJ(

;%
FRM=subs(F,{q1,q2,q3},{Qtem(1),Qtem(2),Qtem(3)});
Qtem=QJ(:,j)-1/2*QQJ(

;
      if(j==2)
         Qtem=QJ(:,j-1);
      end
      FLM=subs(F,{q1,q2,q3},{Qtem(1),Qtem(2),Qtem(3)});
      QJnmid=QJ(:,j)-1/2*tstep/xstep*(FRM-FLM)+1/2*tstep*SJ(:,j);
      QL=QJnmid+1/2*QQJ(

;
Qtem=QJ(:,j+1)+1/2*QQJ1(

;% QQJ1未赋值就用了
      if(j==m)
         Qtem=QJ(:,j+2);
      end
      FRM=subs(F,{q1,q2,q3},{Qtem(1),Qtem(2),Qtem(3)});
      Qtem=QJ(:,j+1)-1/2*QQJ1(

;
      FLM=subs(F,{q1,q2,q3},{Qtem(1),Qtem(2),Qtem(3)});
      QJ1nmid=QJ(:,j+1)-1/2*tstep/xstep*(FRM-FLM)+1/2*tstep*SJ(:,j+1);
      QR= QJ1nmid-1/2*QQJ1(

;%%QQJ1应该可以换成QQJ试试
      lamd_max=max(abs(lamdaJ(:,j)))*max(abs(QR+QL)/2);%需检查
      OLR=tstep/xstep*lamd_max*(QR-QL);
      FR=subs(F,{q1,q2,q3},{QR(1),QR(2),QR(3)});
      FL=subs(F,{q1,q2,q3},{QL(1),QL(2),QL(3)});
      F_TVDLFJ1=1/2*(FL+FR)-1/2*1/(tstep/xstep)*OLR;
      QJn1=QJ(:,j)-tstep/xstep*(F_TVDLFJ1- F_TVDLFJ)+tstep*SJ(:,j);
      F_TVDLFJ=F_TVDLFJ1;
      qq1=QJn1(1);
      qq2=QJn1(2);
      qq3=QJn1(3);
      vEln1(j)=qq2/denL;
      vdenGn1(j)=qq1/(1-vEln1(j));
      vUmn1(j)=qq3/(qq1+qq2);
      if(lamd_max>lamdn_max)
         lamdn_max=lamd_max;
      end
end
vtn=vtn+tstep;
tstep=1*xstep/lamdn_max;
vdenGn1(m+2)=10*100000/(z*R*T);%已知边界条件
vEln1(m+2)=vEln(m+2);
vUmn1(m+2)=vUmn(m+2);
F_TVDLFJ1=subs(F,{q1,q2,q3},{vdenGn1(m+2)*(1-vEln1(m+2)),denL*vEln1(m+2),(vdenGn1(m+2)*(1-vEln1(m+2))+denL*vEln1(m+2))*vUmn1(m+2)});
QJn1=QJ(:,m+1)-tstep/xstep*(F_TVDLFJ1- F_TVDLFJ)+tstep*SJ(:,m+1);
qq1=QJn1(1);
qq2=QJn1(2);
qq3=QJn1(3);
vEln1(m+1)=qq2/denL;
vdenGn1(m+1)=qq1/(1-vEln1(m+1));
vUmn1(m+1)=qq3/(qq1+qq2);
Ptem1=z*R*T*vdenGn(1);
Ptem2=z*R*T*vdenGn(2);
Ptem3=z*R*T*vdenGn1(2);
Ptem=(Ptem3+(Ptem1-Ptem2)+Ptem1)/2;
vdenGn1(1)=Ptem/(z*R*T);
vEln1(m+2)=(vEln1(m+1)+vEln(m+2)-vEln(m+1)+vEln(m+2))/2;
vUmn1(m+2)=(vUmn1(m+1)+vUmn(m+2)-vUmn(m+1)+vUmn(m+2))/2;
n=n+1;
vdenGn=vdenGn1;
vEln=vEln1;
vUmn=vUmn1;
if(mod(n,100)==0)
      clear maplemex;
end
end

[ Last edited by hotoil on 2008-11-20 at 10:12 ]

回复此楼

» 猜你喜欢

1楼 2008-11-19 10:17:32

已阅关注TA 给TA发消息送TA红花 TA的回帖

chjw2010

铁杆木虫 (著名写手)

项目经理

应助: 2 (幼儿园)
金币: 5256.1
散金: 387
红花: 1
帖子: 1517
在线: 97小时
虫号: 522402
注册: 2008-03-10
性别: GG
专业: 有机合成

不清楚，帮楼主顶一下了！

回复此楼

新的一年!

2楼2008-11-19 15:11:51

已阅关注TA 给TA发消息送TA红花 TA的回帖

相关版块跳转我要订阅楼主 hotoil 的主题更新

返回列表