Message Boards Message Boards

Problem using WhenEvent to constrain solution

Julian Lovlie
Posted 12 years ago
Note: This question has also been posted at http://mathematica.stackexchange.com/questions/37907/problem-using-whenevent-to-constrain-solution

Hi guys.

Problem: Simulate pressure in volumes for 1 second. 

The circuit is as follows:












From this I set up the governing DE for both volumes and plot the solution:
 Ad1 = 10/1000^2;
 
 Ad2 = 1.5*1000^(-2);
 
 Ad3 = 1.5*1000^(-2);
 
 Cd1 = 0.67;
 
 Cd2 = 0.67;

Cd3 = 0.67;

V1 = 10/1000;

V2 = 10/1000;

Rho = 875;

beta = 1000*10^6;

ps = 100*10^5;

Q1 = Ad1*Cd1*Sqrt[(2/Rho)*(ps - p1[t])];

Q2 = Ad2*Cd2*Sqrt[(2/Rho)*(ps - p2[t])];

Q5 = Ad3*Cd3*Sqrt[(2/Rho)*(p1[t] - p2[t])];

Q3 = Q1 - Q5;

Q4 = Q2 + Q5;

s = NDSolve[{p1'[t] == (beta*Q3)/V1, p2'[t] == (beta*Q4)/V2,

    p1[0] == p2[0] == 0, WhenEvent[p1[t] >= ps, p1[t] -> ps],

    WhenEvent[p2[t] >= ps, p2[t] -> ps]}, {p1, p2}, {t, 0, 1}];

Plot[Evaluate[{p1[t]/10^5, p2[t]/10^5} /. s], {t, 0, 1}]
Output:




























Mathematica aborts integration after t=0.69 since it encounters complex solutions. This is due to p1 and p2 at some point getting larger than ps, which makes no sense, particularly since I have added 2 'WhenEvent's which shouldn't 'allow' p1 and p2 to be greater than ps (see code). Complex solutions can be avoided by adding 'Abs', however, then the solution seem to completely diverge:



























Question: Why doesn't the 'WhenEvent's seem to 'work'?

PS.; I have obtained a more credible solution using Matlab and the same constrictions:





























Matlab code:

  close all;
  [b][/b]
  clear;
  
  %Basic data
  
  EndTime=1;
  
  StepTime=1e-5;
 
 ps=100*1e5;
 
 Cd1=0.67;
 
 Cd2=0.67;
 
 Cd3=0.67;
 
 Ad1=10*1000^(-2);
 
 Ad2=1.5*1000^(-2);
 
 Ad3=1.5*1000^(-2);
 
 V1=10/1000;
 
 V2=10/1000;
 
 rho=875;
 
 beta=1000*1e6;
 
 
 
 %Initialize
 
 p1_initial=0;
 
 p2_initial=0;
 
 %Initially old values are simply set to current values
 
 Time=0.0;
 
 p1=p1_initial;
 
 p2=p2_initial;
 
 
 
 
 
 %Initialize counters so that plot data is only saved once pr. a number of
 
 %time steps corresponding to ReportInterval
 
 ReportCounter=0;
 
 ReportInterval=10;
 
 Counter=ReportInterval;
 
 %Start time integration
 
 while Time<EndTime
 
 if p1>=ps
 
     p1=ps;
 
 end
 
 if p2>=ps
 
     p2=ps;
 
 end
 
         Q1=Cd1*Ad1*sqrt(2/rho*(ps-p1));
 
         Q2=Cd2*Ad2*sqrt(2/rho*(ps-p2));
 
         Q5=Cd3*Ad3*sqrt(2/rho*(p1-p2));%Flow through orifice 3
 
         Q3=Q1-Q5;%Q3 to volume 1
 
         Q4=Q2+Q5;%Q4 to volume 2
 
         p1Dot=beta*Q3/V1;
 
         p2Dot=beta*Q4/V2;
 
        
 
        
 
         %report
 
         if Counter==ReportInterval

            Counter=0;

            ReportCounter=ReportCounter+1;

            Time_Plot(ReportCounter)=Time;

            p1_Plot(ReportCounter)=p1*1e-5;

            p2_Plot(ReportCounter)=p2*1e-5;

            

            

        end;

        %Time integrate

        p1=p1+p1Dot*StepTime;

        p2=p2+p2Dot*StepTime;

        Time=Time+StepTime;

        Counter=Counter+1;

    end;

    plot(Time_Plot,p1_Plot);

    hold on;

    plot(Time_Plot,p2_Plot,'r');

    grid;
Matlab solution seems to be a good fit to simulated results (from SimulationX):
POSTED BY: Julian Lovlie
3 Replies
Sort By:
I typically would use Re instead of Abs to ensure that Sqrt only gave real values. Sqrt will be problematic if it is given something that gets close to 0 since floating point error could possibly push it below 0 and introduce complex values. A quick solution to this is to use Re with Sqrt:
Q1 = Ad1*Cd1*Re@Sqrt[(2/Rho)*(ps - p1[t])];


Q2 = Ad2*Cd2*Re@Sqrt[(2/Rho)*(ps - p2[t])];


Q5 = Ad3*Cd3*Re@Sqrt[(2/Rho)*(p1[t] - p2[t])];

This seems to give the answer you were looking for.  I'm not familiar with the behavior of WhenEvents in this case, but I would assume it has to do with the fact that the condition would be triggered once p1 had alread crossed ps and not before it. 
POSTED BY: Sean Clarke
Johan Rhodin
Johan Rhodin
Posted 12 years ago
Here is another way to formulate the problem, suggested by Malte:
 Ad1 = 10/1000^2;
 Ad2 = 1.5*1000^(-2);
 Ad3 = 1.5*1000^(-2);
 Cd1 = 0.67;
 Cd2 = 0.67;
 Cd3 = 0.67;
 V1 = 10/1000;
 V2 = 10/1000;
 Rho = 875;
beta = 1000*10^6;
ps = 100*10^5;
Q1 = Ad1*Cd1*Sqrt[Max[(2/Rho)*(ps - p1[t]), 0]];
Q2 = Ad2*Cd2*Sqrt[Max[(2/Rho)*(ps - p2[t]), 0]];
Q5 = Ad3*Cd3*Sqrt[Max[(2/Rho)*(p1[t] - p2[t]), 0]];
Q3 = Q1 - Q5;
Q4 = Q2 + Q5;
eqs = {p1'[t] == (beta*Q3)/V1, p2'[t] == (beta*Q4)/V2,
   p1[0] == p2[0] == 0};
That way you won't run into problems solving it in Mathematica.

If you have a copy of SystemModeler, an equivalent Modelica model would be:
  model complex
   Real p1;
   Real p2;
 initial equation
   p1 = p2;
   p1 = 0;
 equation
   der(p1) = 100000000000 * (0.000006700000000000001 * sqrt(max(0, 2 / 875 * (10000000 + (-1) * p1))) + (-0.0000010050000000000001) * sqrt(max(0, 2 / 875 * (p1 + (-1) * p2))));
   der(p2) = 100000000000 * (0.0000010050000000000001 * sqrt(max(0, 2 / 875 * (10000000 + (-1) * p2))) + 0.0000010050000000000001 * sqrt(max(0, 2 / 875 * (p1 + (-1) * p2))));
end complex;
POSTED BY: Johan Rhodin
Nicolas Venuti
Nicolas Venuti
Posted 12 years ago
Hello Julian

This is a common issue in hydraulic simulation.
The Sqrt function has a vertical tangent at zero,
this make the system undamped which is not realistic,
and drives to a stiff and oscillant solution when approaching null flow.

The realistic way to solve it, is to use a modified sqrt root function (sqr respectively)
This function must integrates two realistic effects :
- the flow can go both ways, so you will not have complex solutions
- the flow does a transition from turbulent to laminar around zero, preserving damping around null flow
ModifiedSignedSqrt[x_] := If[Abs[x]< 1,x,Sign[x]*Sqrt[Abs[x]]]
Plot[ModifiedSignedSqrt[x], {x, -10, 10}]

You must adjust the treshold (here 1) depending of the Reynolds number.
And obviously keep the function continuous around the transition.

Doing this way ensure that the computation is done in reals, and the system is not stiff around null flow.
This will remove the solver issue (event handling is not required), and give a realistic result 
Q1 = Ad1*Cd1*ModifiedSignedSqrt[(2/Rho)*(ps - p1[t])]
...
eqs = {p1'[t] == (beta*Q3)/V1, p2'[t] == (beta*Q4)/V2,    p1[0] == p2[0] == 0};
...
PS : If you are interested by realistic modeling, try bond graphs.
I am the author of the bond graph toolbox for mathematica

Nicolas
POSTED BY: Nicolas Venuti
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract