Note: This question has also been posted at the Wolfram Community
Problem: Simulate pressure in volume 1 and 2 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}]
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 WhenEvent
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 an OK fit with simulation results (from SimulationX):
Subscript
occurrences and do not use:=
to set values to simple constants use=
. By the wayC
is reserved for something special in Mathematica, check the doc. $\endgroup$