Search All of the Math Forum:
Views expressed in these public forums are not endorsed by
NCTM or The Math Forum.


Torsten
Posts:
1,681
Registered:
11/8/10


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


"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.649647e01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.776357e15) 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 * ((1t)/design.t1); > else > Av = 0; > end > > % Working out the flow rate through either side of the valve. > Q(1,1) = sign(design.poz(1)) * design.cd * Av * ((2*abs(design.poz(1))/design.rho)^(1/2)); > > Q(1,2) = sign(design.psz(2)) * design.cd * Av * ((2*abs(design.psz(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.7e6; % Volume from the valve to the cylinder (in pipework) > design.la = 0.5; % Stroke length of the cylinder > design.lsl = 0.01; > design.lgap = 20e6; % Anulus gap > design.w = 0.01; % Valve width opening > design.xvmax = 0.5e3; % 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 = 15e3; % 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 0t1 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 1e6. > %odeoptions = odeset('RelTol', 1e3); > 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 ifstatement 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 * ((1t)/design.t1); else Av = 0; end 5.The sign and absfunctions may introduce problems because of their nondifferentiability.
Good luck!
Best wishes Torsten.



