Message Boards Message Boards

WhenEvent stops working in the middle of solution of ODEs

Posted 8 years ago

I am solving a set of ODEs and equations to determine temporal variation of Temperature, Pressure and Volume in a real gas, with discrete decision points based on WhenEvent. The solution works fine at start. However the WhenEvent logic stops working all of a sudden after time step t=78.26 (see plot of Pressure vs time below). I would appreciate any help or suggestion to identify the cause of WhenEvent breaking down.

enter image description here

Clear["Global`*"];
dens = ThermodynamicData["Nitrogen", 
   "Density", {"Pressure" -> Quantity[17, "MegaPascals"], 
    "Temperature" -> Quantity[300, "Kelvins"]}];
\[Rho] = 184.;
M = 28.013; v = M/\[Rho];
tend = 100;
f = .1; tp = 1/f; \[Omega] = 2 \[Pi] f; mg = V0*\[Rho]; Cv = .743;
D0 = 250 10^-3; L = 4.5;
V0 = \[Pi]/4 D0^2 L;
Aw = \[Pi] D0 L ;
h = .575;
\[Tau] = (mg Cv)/(h Aw) ;
Tw = 300;
R = 8.314;
a = 2.54; A0 = 106.73; b = .002328; B0 = .04074; c = 7.379 10^4; C0 = 
 8.164 10^5; \[Alpha] = 1.272 10^-4; \[Gamma] = .0053;

ul = 100;
vc = .005;
pmin = 10 10^6;
pmax = 30 10^6;
Pwr = 50 10^3;

sol = Quiet@
   NDSolve[{Temp'[t] - (Tw - Temp[t])/\[Tau] + 
       1/Cv ((R Temp[t])/v (1 + b/v^2) + 
          1/v^2 (B0 R Temp[t] + 2 C0/(Temp[t])^2) - (2 c)/(
           v^3 (Temp[t])^2) (1 + \[Gamma]/v^2) Exp[-(\[Gamma]/v^2)])*
        Vol'[t] == 0, Press[t]*Vol'[t] + ax[t]*Pwr == 0, 
     Press[t] == ((R Temp[t])/v + 
         1/v^2 (B0 R Temp[t] - A0 - C0/(Temp[t])^2) + 
         1/v^3 (b R Temp[t] - a) + 1/v^6 a \[Alpha] + 
         1/(v^3 (Temp[
             t])^2) (c (1 + \[Gamma]/v^2) Exp[-(\[Gamma]/v^2)])) * 
       10^3 , Temp[0] == 300, Vol[0] == V0, Press[0] == pmin, 
     ax[0] == 1, 
     WhenEvent[
      Press[t] > pmax && Press'[t] > 0, {t1 = t, Print["t1 = ", t1], 
       ax[t] -> 0}], 
     WhenEvent[
      t > t1 + 10, {t2 = t, Print["t2 = ", t2], ax[t] -> -1}], 
     WhenEvent[
      Vol[t] == V0  && Press'[t] < 0 , {t3 = t, Print["t3 = ", t3], 
       ax[t] -> 0}], 
     WhenEvent[
      t > t3 + 10, {t4 = t, Print["t4 = ", t4], ax[t] -> 1}]}, {Temp, 
     Vol, Press}, {t, 0., tend}, DiscreteVariables -> ax, 
    AccuracyGoal -> 8, PrecisionGoal -> 8];
Plot[Evaluate[Press[t] /. sol] 10^-6, {t, 0, tend}, Frame -> True, 
 Axes -> True, PlotRange -> All, 
 FrameLabel -> {"Time (secs)", "Pressure (MPa)"}, 
 GridLines -> Automatic, AspectRatio -> .65, 
 PlotStyle -> {{Brown, Thickness[0.0075]}}, 
 BaseStyle -> {FontWeight -> "Bold", FontSize -> 15}]

Plot[Evaluate[Temp[t] /. sol], {t, 0, tend}, Frame -> True, 
 Axes -> True, PlotRange -> All, 
 FrameLabel -> {"Time (secs)", "Temperature (K)"}, 
 GridLines -> Automatic, AspectRatio -> .65, 
 PlotStyle -> {{Blue, Thickness[0.01]}}, 
 BaseStyle -> {FontWeight -> "Bold", FontSize -> 15}]

Plot[Evaluate[Vol[t] /. sol], {t, 0, tend}, Frame -> True, 
 Axes -> True, PlotRange -> All, 
 FrameLabel -> {"Time (secs)", "Volume (m3)"}, GridLines -> Automatic,
  AspectRatio -> .65, PlotStyle -> {Gray, Thickness[0.01]}, 
 BaseStyle -> {FontWeight -> "Bold", FontSize -> 15}]
POSTED BY: Kaushik Mallick

Can anyone help me please?

POSTED BY: Kaushik Mallick
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract