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