#*-------------------------------------------------------------------
* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
*
* This LIBRARY is free software; you can distribute it and/or modify
* it under the therms of the ALSOC FREE LICENSE as available at
* http://www.enq.ufrgs.br/alsoc.
*
* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
* from http://www.rps.eng.br Copyright (C) 2002-2004.
* All rights reserved.
*
* EMSO is distributed under the therms of the ALSOC LICENSE as
* available at http://www.enq.ufrgs.br/alsoc.
*
*--------------------------------------------------------------------
* The car axis problem
*
* The problem is a stiff DAE of index 3, consisting of 8
* differential and 2 algebraic equations.
*--------------------------------------------------------------------
* Author: Rafael de Pelegrini Soares
* $Id: sample_car.mso 83 2006-12-08 20:29:34Z paula $
*--------------------------------------------------------------------*#
FlowSheet CarAxis
PARAMETERS
l as Real(Default=1);
l0 as Real(Default=1/2);
eps as Real(Default=1e-2);
M as Real(Default=10);
h as Real(Default=1/5);
tau as Real(Default=3.1514/5);
w as Real(Default=10);
r as Real(Default=0.1);
VARIABLES
xl as Real(Default = 0, Lower=-1, Upper = 1);
yl as Real(Default = 0.5, Lower=-1, Upper = 1);
xr as Real(Default = 1, Lower=-2, Upper = 2);
yr as Real(Default = 0.5, Lower=-1, Upper = 1);
q(4) as Real(Default = -1);
lambda(2) as Real(Default = 0);
xb as Real(Default = 1);
yb as Real(Default = 0);
ll as Real(Default = 0.5, Brief="Left spring length");
lr as Real(Default = 0.5, Brief="Right spring length");
EQUATIONS
diff([xl, yl, xr, yr]) = q;
eps^2*M/2 * diff(q(1)) = (l0-ll)*xl/ll
+lambda(1)*xb +2*lambda(1)*(xl-xr);
eps^2*M/2 * diff(q(2)) = (l0-ll)*yl/ll
+lambda(1)*yb +2*lambda(2)*(yr-yl) - eps^2*M/2;
eps^2*M/2 * diff(q(3)) = (l0-lr)*(xr-xb)/lr
-2*lambda(2)*(xl-xr);
eps^2*M/2 * diff(q(4)) = (l0-lr)*(yr-yb)/lr
-2*lambda(2)*(yl-yr) - eps^2*M/2;
xl*xb + yl*yb = 0;
(xl-xr)^2 + (yl-yr)^2 = l^2;
xb = sqrt(l^2-yb^2);
yb = r*sin(w*time);
ll = sqrt(xl^2 + yl^2);
lr = sqrt((xr-xb)^2 + (yr-yb)^2);
INITIAL
yl = 0.5;
xr = 1;
q(3) = -0.5;
q(4) = 0;
OPTIONS
time = [0:0.01:3];
integration = "index0"; # "original";
#DAESolver = "mebdf";
relativeAccuracy = 1e-5;
absoluteAccuracy = 1e-5;
indVarAccuracy = 1e-3;
end