Path: news.mathworks.com!not-for-mail
From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: ODE stiffness problems
Date: Thu, 14 Mar 2013 08:52:06 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 129
Message-ID: <khs33m$h41$1@newscl01ah.mathworks.com>
References: <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 1363251126 17537 172.30.248.37 (14 Mar 2013 08:52:06 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Thu, 14 Mar 2013 08:52:06 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 3799640
Xref: news.mathworks.com comp.soft-sys.matlab:791119

"Matthew" wrote in message <khqcg9$anl$1@newscl01ah.mathworks.com>...
> 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

1. Call the integrator as
[t, z] = ode15s(@(t,z)cylinderoderhs(t,z,design), span, [sp1 sp2 sxa sva], odeoptions);
2. Define your function as
 function dz = cylinderoderhs(t,z,design)
3.Insert Q(1,1) and Q(1,2) in the calculation of dz(1) and dz(2) instead of Q(1) and 
Q(2).
4. For clarity, change your if-statement in the calculation of Av to
 if t<design.t1 
     Av = design.w * design.xvmax * (t/design.t1);
 elseif t>= design.t1 & t<design.t2
     Av = design.w * design.xvmax;
 elseif t>= t.design.t2 & t<design.t3
     Av = design.w * design.xvmax * ((1-t)/design.t1);
 else
     Av = 0;
 end
5.The sign and abs-functions may introduce problems because of their non-differentiability.

Good luck!

Best wishes
Torsten.