Math Forum Discussions

Torsten

Posts: 1,136
Registered: 11/8/10

Re: ODE stiffness problems
Posted: Mar 14, 2013 4:52 AM

Plain Text

"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.