1
|
20774 Views
|
3 Replies
|
4 Total Likes
View groups...
Share
GROUPS:

# Problem using WhenEvent to constrain solution

Posted 11 years ago
3 Replies
Sort By:
Posted 11 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[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 mathematicaNicolas
Posted 11 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 11 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.