MATLAB Answers

Reduce computing time ode system

조회 수: 4(최근 30일)
EldaEbrithil
EldaEbrithil 21 Jun 2020
편집: EldaEbrithil 22 Jun 2020
Hi all
i would like to reduce the computing time of my code in such a way as to increase the dimensions of integration domain (L), hradialchmaber vetor, and M01 vector
hradialchamber=(1.5:0.5:2);%mm
for i=1:length(hradialchamber)
Hradialchamber(i)=hradialchamber(i)/1000;%m
end
rextbobbin=9.525:0.3:17.5;%
Rextbobbin=rextbobbin/1000;%m
Atomicweight=131.29;%
gamma=1.667;%
R=8314/Atomicweight;%
Cp=(gamma/(gamma-1))*R;%e
dconduttore=(0.43);%mm
Dconduttore=dconduttore/1000;%m
throatRadius=0.100:0.0005:0.5;%mm
throatradius=throatRadius/1000;
exitRadius=0.600:0.001:3;%mm
exitradius=exitRadius/1000;
exitarea=pi*exitradius.^2;%m^2
Mexit=5.88;%INPUT
Tcoldgas=300;%K
for i=1:length(throatRadius)
for j=1:length(exitRadius)
radiusparameter=0.5;
Mcheck(i,j)=((0.50*exitradius(j)/(throatradius(i)))^2)-...
((1/Mexit)*((2/(gamma+1))*...
(1+(((gamma-1)/2)*Mexit^2)))^((gamma+1)/(2*(gamma-1))));%
diffM(i,j)=Mcheck(i,j)-Mexit;
end
end
DiffminM2=min(diffM(diffM>0));
diffminM2=max(diffM(diffM<0));%
if DiffminM2-abs(diffminM2)>0%
[righ2,colh2]=find(diffM==diffminM2);
for i=1:length(righ2)
Mcheck1(i)=Mcheck(righ2(i),colh2(i));
end
exitradius1=exitradius(colh2);
throatradius1=throatradius(righ2);
else DiffminM2-abs(diffminM2)<0;
[RigH,ColH]=find(diffM==DiffminM2);
for i=1:length(RigH)
Mcheck1(i)=Mcheck(RigH(i),ColH(i));
end
exitradius1=exitradius(ColH);
throatradius1=throatradius(RigH);
end
exitradiusfinal=max(exitradius1);
throatradiusfinal=max(throatradius1);
throatarea=pi*throatradiusfinal^2;%m^2
exitareafinal=pi*exitradiusfinal^2;%m^2
pressure=2.2;%bar
Pressure=pressure*10^5;%Pa
theater=803;%K
emittancetantalum=0.16;
dexttantalum=1.5;%mm
Dexttantalum=dexttantalum/1000;%m
for i=1:length(hradialchamber)
Aeffettiva(i)=(Dexttantalum*(Hradialchamber(i)))-(pi*(Dexttantalum^2)/4);
A2(i)=Aeffettiva(i);
M2=0.001:0.0001:0.8;
for j=1:length(M2)
M2check(j,i)=(A2(i)/(throatarea))-...
((1/M2(j))*((2/(gamma+1))*...
(1+(((gamma-1)/2)*M2(j)^2)))^((gamma+1)/(2*(gamma-1))));
end
end
Diff_minM2=min(M2check(M2check>0));
diff_minM2=max(M2check(M2check<0));
if Diff_minM2-abs(diff_minM2)>0%
[righ_2,colh_2]=find(M2check==diff_minM2);
M2final=M2(righ_2);
else Diff_minM2-abs(diff_minM2)<0
[RigH_2,ColH_2]=find(M2check==Diff_minM2);
M2final=M2(RigH_2);
end
L=0.2:0.1:0.3;
M01=(0.1:0.1:M2final);
dynamicviscosity=54.89*10^-6;
dynamicviscositycold=23.23*10^-6;
dynamicviscosityaverage=(dynamicviscosity+dynamicviscositycold)/2;
thermalconductivityxenonhot=13.07*10^-3;
thermalconductivityxenoncold=5.555*10^-3;
thermalconductivityxenonaverage=(thermalconductivityxenonhot+...
thermalconductivityxenoncold)/2;
Pr=dynamicviscosityaverage*Cp/thermalconductivityxenonaverage;
ro1=Pressure/(R*Tcoldgas);
for y=1:length(hradialchamber)
Aeffettiva(y)=(Dexttantalum*(Hradialchamber(y)))-(pi*(Dexttantalum^2)/4);
Lq(y)=(Aeffettiva(y)^0.5);
Perimetro(y)=Hradialchamber(y)*2+Dexttantalum*2+(2*pi*(Dexttantalum/2));
Dh(y)=(4*(Aeffettiva(y))/Perimetro(y));
sqrtAeffettiva(y)=Aeffettiva(y)^0.5;
for k=1:length(M01)
u1(k)=M01(k)*(gamma*R*Tcoldgas)^0.5;
mdotstart(k,y)=ro1*u1(k)*Aeffettiva(y);
Re(k,y)=(4*mdotstart(k,y))/(dynamicviscosityaverage*pi*Dh(y));
Deannumber(k,y)=Re(k,y)*(thermalconductivityxenonaverage)^0.5;
for i=1:length(rextbobbin)
delta(i)= Dexttantalum/(Rextbobbin(i));
Recritic(i)=2100*(1+12*(delta(i)^0.5));
Rediffer(i,k,y)=Re(k,y)-Recritic(i);
%%%%From here to end the critical part%%%
if Re(k,y)>Recritic(i)
RE(k,y)=Re(k,y);
[RErig,REcol]=find(RE>0);
A(y)=Aeffettiva(y);
Mo1(k)=M01(k);
NuT(k,y)=Re(k,y)*(Pr^(1/3))*(1/(2*(3.6*log10(6.115/Re(k,y)))^2));
fT(k,y)=0.0767/(Re(k,y)^0.25);
hT(k,y)=NuT(k,y)*thermalconductivityxenonaverage/(sqrtAeffettiva(y));
ChT(k,y)=NuT(k,y)/(Re(k,y)*Pr);
end
end
end
end
Nc=0.001;
P0=2.200000e+05;
T0=300;
Tt0=T0;
Pt0=P0;
for y=1:length(hradialchamber)
Aeffettiva_cell(y)={Aeffettiva(y)};
Perimetro_cell(y)={Perimetro(y)};
for j=1:length(M01)
Y0(j)={[Tt0,Pt0,Mo1(j)]};
ChT_cell(j,y)={ChT(j,y)};
ft_cell(j,y)={fT(j,y)};
end
end
for i=1:length(L)
Dall(i) ={[0:Nc:L(i)]};
end
for iDom = 1:numel(Dall)
xRange = Dall{iDom};
for iInitial = 1:numel(Y0)
for iAeff=1:numel(Aeffettiva_cell)
Ch=ChT_cell{iInitial,iAeff};%it doesn't depend by iDom
Fc=ft_cell{iInitial,iAeff};
Aeff=Aeffettiva_cell{iAeff};
Perim=Perimetro_cell{iAeff};
[xSol{iDom,iInitial,iAeff},YSol{iDom,iInitial,iAeff}]=ode23(@(x,Y) ...
chambsinglebobb(x,Y,Ch,Aeff,Perim,Fc) ,xRange,Y0{iInitial});
end
end
end
function dYdx=chambsinglebobb(x,Y,Ch,Aeff,Perim,Fc)
gamma=1.667;
Dexttantalum=0.001500000000000;
Tw=803;
Tt=Y(1);
Pt=Y(2);
M=Y(3);
dTtdx=Ch*(Tw-Tt)*(Perim/Aeff);
dPtdx=Pt*((-gamma*((M^2)/2)*Fc*(Perim/Aeff))-...
(gamma*((M^2)/2)*dTtdx*(1/Tt)));
dMdx=M*(((1+((gamma-1)/2)*M^2)/(1-M^2))*((gamma*(M^2)*Fc*...
Perim)/(2*Aeff))+...
((0.5*((gamma*M^2)+1))*(1+((gamma-1)/2)*M^2)/(1-M^2))*...
(Ch*(Tw-Tt)*Perim)/(Aeff*Tt));
dYdx=[dTtdx;dPtdx;dMdx];
end
In the order, the operations of the code are the calculation of:
  • throatradiusfinal
  • exitradiusfinal
  • throatarea
  • M2final
  • NuT
  • ChT
  • fT
  • YSol
The variable parameters are: hardialchamber, rextbobbin, M01, L. The problem related to high computational time is connected (i think) with the the calculation of the ode system solution. In particular i have written the code in such a way the ODE system is solved numerous time namely for different IC (showbelow), different value of Ch in the differential equations, different value of Fc in the differential equations and of course different integration domain (L):
Y0(j)={[Tt0,Pt0,Mo1(j)]}
Ch=ChT_cell{iInitial,iAeff};
Fc=ft_cell{iInitial,iAeff};
xRange = Dall{iDom};
The problem (i think) is connected with ChT that it is physically independent of L, this leads to increase very much the number of combinations allowed. I have written the same code but with the only difference of ChT dependent on L, this code is much faster than the privious one BUT of course it is not physically acceptable. So i ask you if you know some shrewdness for increasing the calculation time for this problem.
Regards

답변(2개)

Alan Stevens
Alan Stevens 21 Jun 2020
Here's a start, looking at your triply nested y, k and i loops.
% Remove expressions that don't depend on y, k or i from the loops
dynamicviscosity=54.89*10^-6;
dynamicviscositycold=23.23*10^-6;
dynamicviscosityaverage=(dynamicviscosity+dynamicviscositycold)/2;
thermalconductivityxenonhot=13.07*10^-3;
thermalconductivityxenoncold=5.555*10^-3;
thermalconductivityxenonaverage=(thermalconductivityxenonhot+...
thermalconductivityxenoncold)/2;
Pr=dynamicviscosityaverage*Cp/thermalconductivityxenonaverage;
ro1=Pressure/(R*Tcoldgas);
for y=1:length(hradialchamber)
% Put expressions that depend only on y here
Aeffettiva(y)=(Dexttantalum*(Hradialchamber(y)))-(pi*(Dexttantalum^2)/4);
Lq(y)=(Aeffettiva(y)^0.5);
Perimetro(y)=Hradialchamber(y)*2+Dexttantalum*2+(2*pi*(Dexttantalum/2));
Dh(y)=(4*(Aeffettiva(y))/Perimetro(y));
sqrtAeffettiva(y)=Aeffettiva(y)^0.5;
for k=1:length(M01)
% Put expressions that depend on k or k and y here
u1(k)=M01(k)*(gamma*R*Tcoldgas)^0.5;
mdotstart(k,y)=ro1*u1(k)*Aeffettiva(y);
Re(k,y)=(4*mdotstart(k,y))/(dynamicviscosityaverage*pi*Dh(y));
Deannumber(k,y)=Re(k,y)*(thermalconductivityxenonaverage)^0.5;
for i=1:length(rextbobbin)
% Put expressions that depend on i, i and y, i and k or i, y and k here
delta(i)= Dexttantalum/(Rextbobbin(i));
Recritic(i)=2100*(1+12*(delta(i)^0.5));
Rediffer(i,k,y)=Re(k,y)-Recritic(i);
%%%%From here to end the critical part%%%
if Re(k,y)>Recritic(i)
RE(k,y)=Re(k,y);
[RErig,REcol]=find(RE>0);
A(y)=Aeffettiva(y);
Mo1(k)=M01(k);
NuT(k,y)=Re(k,y)*(Pr^(1/3))*(1/(2*(3.6*log10(6.115/Re(k,y)))^2));
fT(k,y)=0.0767/(Re(k,y)^0.25);
hT(k,y)=NuT(k,y)*thermalconductivityxenonaverage/(sqrtAeffettiva(y));
ChT(k,y)=NuT(k,y)/(Re(k,y)*Pr);
end
end
end
end
  댓글 수: 1
EldaEbrithil
EldaEbrithil 21 Jun 2020
Thank you for reply me
ok i have done it, but the computational time is still too high
if you have any question please let me know

댓글을 달려면 로그인하십시오.


Alan Stevens
Alan Stevens 21 Jun 2020
Look at the other triple loop. Similar improvements can be made:
for y=1:length(hradialchamber)
Aeffettiva_cell(y)={Aeffettiva(y)};
Perimetro_cell(y)={Perimetro(y)};
for j=1:length(M01)
Y0(j)={[Tt0,Pt0,Mo1(j)]};
ChT_cell(j,y)={ChT(j,y)};
ft_cell(j,y)={fT(j,y)};
end
end
for i=1:length(L)
Dall(i) ={[0:Nc:L(i)]};
end
  댓글 수: 6
EldaEbrithil
EldaEbrithil 22 Jun 2020
ok i have found the problem: the ode23 solve the system which gives some values of
YSol{iDom,iInitial}(iDom,3)
bigger or equal to 1. This is not allowed by the equations of the system, this create a lot of problems in the solving process and therefore the long computational times are explained. Is there a way to prevent ode23 to solve the system as soon as
YSol{iDom,iInitial}(iDom,3)
values are near 1 or maybe bigger than one?
something like that but i don't know:
if YSol{i,j}(i,3)>1
break;
end

댓글을 달려면 로그인하십시오.

제품


릴리스

R2019a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by