# Problem using WhenEvent to constrain solution

Posted 10 years ago
19778 Views
|
3 Replies
|
4 Total Likes
|
 Note: This question has also been posted at http://mathematica.stackexchange.com/questions/37907/problem-using-whenevent-to-constrain-solutionHi 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 == p2 == 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=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): 3 Replies
Sort By:
Posted 10 years ago
 Hello JulianThis 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 flowModifiedSignedSqrt[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 == p2 == 0};...PS : If you are interested by realistic modeling, try bond graphs.I am the author of the bond graph toolbox for mathematicaNicolas
Posted 10 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 == p2 == 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 10 years ago
 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.