Hi all
I have solved a system of differential equations in an iterative way, as you can see from the graph I have several curves each of which is associated with a different combination of the iterative parameters. Again from the graph, it can be seen that some curves at a certain abscissa diverge; this is not acceptable so I would like to remedy it. Is it possible to create a new solution vector capable of avoiding the divergence as described below ?
I would like the solution not to diverge near the various cusp points and, on the contrary, to maintain the value immediately upstream of the latter, until the end of the domain. I'll explain: suppose we have the following solution values 1000 K,1020 K,2000 K,4000 K I would like instead: 1000 K,1020 K,1020 K,1020 K
Thank you for the help
Regards
T0=293;delta_Thot=2500;%s
%%%%%%%%%%%%%
Nc_Thot=2;[ya,sa,lea]=size(Aexposehot);for y=1:yV%passo,altezza tubo,Rextbobbin,lunghezza cavo
for le=1:leV%lunghezza tubo equilibrio
for sv=1:sV%spessore parete3
% for k=1:km%Mach
for rf=1:length(Pressure_percentage_hot)Rextbobb_cellhot(y)={Rextbobbin_fluxhot(y)};Y0hot(y,le)={[T0,m0hot(y,le)]};Dall_Th_hot ={[0:Nc_Thot:delta_Thot]}; Lequilibrium_cellhot(le)={Lequilibriumhot(le)}; Passo_cellhot(y)={Passo_fluxhot(y)}; hrad_cellhot(y)={Hradial_fluxhot(y)};Spessore_parete3_cellhot(sv) ={Spessore_parete3hot(sv)};Aexpose_cellhot(y,sv,le)={Aexposehot(y,sv,le)};Volume_cell_Aisihot(le,y,sv)={Volume_Aisihot(le,y,sv)};Volume_cell_gashot(y,le)={Volume_gashot(y,le)};throat_cell_hot(y)={ThrtR_fluxhot(y)};Pressure_end_chamber_hot_cell(rf)={Pressure_end_chamber_hot(rf)}; endend end endfor iDom_Th = 1:numel(Dall_Th_hot) for leLE=1:numel(Lequilibrium_cellhot) for yY=1:numel(hrad_cellhot) for rF=1:numel(Pressure_end_chamber_hot_cell) for sS=1:numel(Spessore_parete3_cellhot)TSol_hot(iDom_Th,yY,leLE,sS,rF)={zeros(length(Dall_Th_hot{iDom_Th}),2)}; end end end end end opt_t_hot=odeset('NormControl','Refine','Stats','MaxStep');for iDom_Th = 1:numel(Dall_Th_hot) tRange_hot = Dall_Th_hot{iDom_Th}; for leLE=1:numel(Lequilibrium_cellhot) for yY=1:numel(hrad_cellhot) for rF=1:numel(Pressure_end_chamber_hot_cell) for sS=1:numel(Spessore_parete3_cellhot) Ltube_hot= Lequilibrium_cellhot{leLE};%L_cellhot{yY};
Leq_hot= Lequilibrium_cellhot{leLE} ; Vol_Aisi_hot=Volume_cell_Aisihot{leLE,yY,sS}; Aexp_hot=Aexpose_cellhot{yY,sS,leLE}; Hradial_hot=hrad_cellhot{yY}; ThroatradiuS_hot=throat_cell_hot{yY}; Pressure_hot_end_chamber=Pressure_end_chamber_hot_cell{rF}; if Vol_Aisi_hot>0&&Leq_hot>0 &&Aexp_hot>0 &&Ltube_hot>0 && ThroatradiuS_hot>0&&Pressure_hot_end_chamber>0 [tSol{iDom_Th,yY,leLE,sS,rF},TSol_hot{iDom_Th,yY,leLE,sS,rF}]=ode23t(@(t,Y) ... temperaturaheaterclassica(t,Y,Ltube_hot,Voltage1,Vol_Aisi_hot,Aexp_hot,ThroatradiuS_hot,Pressure_hot_end_chamber),tRange_hot,Y0hot{yY,leLE},opt_t_hot); if(any(TSol_hot{iDom_Th,yY,leLE,sS,rF}(:,1)>1600)) % disp([ ' condizione verificata'])
%disp (iDom_Th)
disp(yY) disp(leLE) %disp(sS)
%disp(rF)
break end end end end end endend %end
[tT,yT,leT,sT,rFT]=size(TSol_hot);for t = 1:tT for y=1:yT for le=1:leT for s=1:sT for rf=1:rFT if TSol_hot{t,y,le,s,rf}(end,1)>0 figure(1);set(gcf,'Visible', 'on')plot(tSol{t,y,le,s,rf}(:,1),TSol_hot{t,y,le,s,rf}(:,1))xlabel('time [s]')ylabel('Teq-transition/turbulent [K]') hold onendendendendendend
Best Answer