Path: news.mathworks.com!not-for-mail
From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: ODE stiffness problems
Date: Wed, 13 Mar 2013 17:20:09 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 105
Message-ID: <khqcg9$anl$1@newscl01ah.mathworks.com>
Reply-To: <HIDDEN>
NNTP-Posting-Host: www-05-blr.mathworks.com
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: newscl01ah.mathworks.com 1363195209 10997 172.30.248.37 (13 Mar 2013 17:20:09 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Wed, 13 Mar 2013 17:20:09 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 3714708
Xref: news.mathworks.com comp.soft-sys.matlab:791074

Hi, I am working on a project simulating a moving hydraulic cylinder with fluid being pushed out and drawn in through different ends. The fluid moves through valves which are opening and shutting.

I didn't know anything about functions and ODEs until recently, when someone tried to help me solve the problem. However they do not know what to do now. 

I'm getting this error 'Warning: Failure at t=9.649647e-01.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.776357e-15) at time t. '

There are 4 integrals which need to be calculated by the ODE, position of the piston head, velocity of the piston head, pressure inside side A of the chamber and pressure inside chamber B of the cylinder.

I am following an example that I found, and I get the correct answers up until 1 second. When I simulate more than 1 second I get the error. I have tried messing with the tolerances and using all the different ODE solvers. Can someone please help me. Thanks.

The function is listed below followed by the main script.

function [dz, Av, Q] = cylinderoderhs(t,z,design)

% Calculating the size of the valve open area as a function of time
if t<design.t1
    Av = design.w * design.xvmax * (t/design.t1);
elseif t<design.t2
    Av = design.w * design.xvmax;
elseif t<=design.t3
    Av = design.w * design.xvmax * ((1-t)/design.t1);
else
    Av = 0;
end

% Working out the flow rate through either side of the valve. 
Q(1,1) = sign(design.po-z(1)) * design.cd * Av * ((2*abs(design.po-z(1))/design.rho)^(1/2));

Q(1,2) = sign(design.ps-z(2)) * design.cd * Av * ((2*abs(design.ps-z(2))/design.rho)^(1/2));

dz = zeros(4,1);

% Calculating the change in pressure in chamber side 1, dp1/dt 
dz(1) = (design.b / (design.Vt - design.Acap * z(3))) * (Q(1) + design.Acap * z(4));

% Calculating the change in pressure in chamber side 2, dp2/dt
dz(2) = (design.b /(design.Vt + design.Acap * z(3))) * (Q(2) - design.Acap * z(4));

% Calculating the position of the cylinder, xa. Which is the integrated
% va.
dz(3) = z(4);

% calculating the velocity of the cylinder, va. 
dz(4) = (z(2) * design.Acap - z(1)*design.Arod - design.c*z(4) - design.ma*design.g) / design.ma;

end




% The times for the simulation start and end.
span = [0 1];

% System parameters
design.ma = 750;     % Mass
design.da = 0.05;    % Diameter of the piston
design.drod = 0.025; % Diameter of the rod
design.V = 15.7e-6;  % Volume from the valve to the cylinder (in pipework)
design.la = 0.5;     % Stroke length of the cylinder
design.lsl = 0.01;
design.lgap = 20e-6; % Anulus gap
design.w = 0.01;     % Valve width opening
design.xvmax = 0.5e-3;  % Maximum valve opening (Vertical)
design.cd = 0.62;    % Orafice discharge coefficient
design.g = 9.81;     % Gravity
design.rho = 855;    % Density
design.b = 1.2e9;    % Bulk modulus
design.mu = 15e-3;   % Absolute viscosity
design.ps = 20e6;    %System pressure
design.po = 0;       %Drain pressure

% Annulus areas are the same.
design.Acap = (design.da^2 * pi/4) - (design.drod^2*pi/4);
design.Arod = design.Acap;

% The total volume in each side of the cylinder at x=0 (Center)
design.Vt = design.V + (design.Acap * design.la)/2;

% The damping coefficent.
design.c = pi * design.da * design.lsl * design.mu / design.lgap;

% Starting postion, xa, velocity, va and pressures either side of the
% chamber.
sxa = 0;
sva = 0;
sp2 = 750*9.81/(0.05^2*pi/4-(0.025^2*pi/4));
sp1 = 0;

% Times that the valve opens. From 0-t1 the valve is opening, it is open
% fully until t2 is reached then takes 0.05 seconds to close.
design.t1 = 0.05;
design.t2 = 0.95;
design.t3 = 1;

% Altering the relative tolerance of the ODE from default 1e-6.
%odeoptions = odeset('RelTol', 1e-3);
odeoptions = [];

% solve the ode function.
[t, z] = ode15s(@cylinderoderhs, span, [sp1 sp2 sxa sva], odeoptions, design);

% Extracting the internally calculated values
for i =1:numel(t)
    [dz(i,:),  Av(i,1), Q(i,1:2)] = cylinderoderhs(t(i), z(i,:), design);
end