Archived OpenModelica forums. Posting is disabled.

Alternative forums include GitHub discussions or StackOverflow (make sure to read the Stack Overflow rules; you need to have well-formed questions)


Forgot password? | Forgot username? | Register

Best practices for modeling rigid stops with Modelica

Best practices for modeling rigid stops with Modelica

Good morning.
I'm trying to model a spool that runs into a valve body.
The stroke of the spool is mechanically limited to the positions x1 = 0 and x2 = 3.3 mm (millimeters).
There is a spring, which pushes the spool to x1.
By an external force Fm applied to the spool, this one should therefore move between x1 <= x <= x2.
I have so a problem of rigid constraints.
I have modeled the rigid constraints with a viscoelastic model (Kelvin-Voigt) and it works quite well. For example with Fm = 1000 N, x is about 0.5 mm (correct).
If Fm = 0 x is around zero (-0.0015, this depends on the elasticity) if Fm = 10000 N (for example) the stroke is stopped at about 3.3 mm.
Here's the model I wrote:

model VS30CC
  parameter Real m(unit="m") = 1;
  parameter Real k(unit="N/mm") = 192.6, L0(unit="mm") = 38;
  parameter Real c(unit="N.s/m") = 50;
  parameter Real d(unit="mm") = 6;
  parameter Real Ptar(unit="bar")=320;
  parameter Modelica.SIunits.Force Fm = 1000;
  Real x(unit="mm",start=0),v(unit="m/s",start=0);
  Modelica.SIunits.Force F,Fk;
  protected
  parameter Real x1(unit="mm")=0,x2(unit="mm")=3.3;
  parameter Real E(unit="N/mm2")=210000;
  parameter Real A(unit="mm2")=0.25*Modelica.Constants.pi*d^2;
  parameter Real L(unit="mm")=10;
  parameter Real Kv(unit="N/mm") = E*A/L;
  parameter Real Ck(unit="N.s/m") = -1000;
  parameter Modelica.SIunits.Force F0 = 0.1*Ptar*A;
equation
  if x<x1 then
    Fk = -Kv*(x-x1)+Ck*v;
  elseif x>x2 then
    Fk = -Kv*(x-x2)+Ck*v;
  else
    Fk = 0;
  end if;
  v = der(0.001*x);
  m * der(v) = F+Fk;
  F = Fm-c*v-k*x-F0;
end VS30CC;

The problem is that I have already tried this method for problems like this, but more complex, where there are forces of surfaces variable in time (the result of pressure variables) and my experience is that the calculation becomes very very slow due the stiffness around the constraints. I submit the basic example and not the more complex because the other codes are very long and not focused to the problem of the constraints. In those codes I also tried to round Fk with the hyperbolic tangent (tanh) but the execution is still very slow.

So I wanted to try something like the classic example "BouncingBall" but with both constraints, so I adjusted the example I found on the Tiller book "Introduction to physical modeling with Modelica". Here's how I've adapted my model:

model VS30CC
  parameter Real m(unit="m") = 1;
  parameter Real k(unit="N/mm") = 192.6, L0(unit="mm") = 38;
  parameter Real c(unit="N.s/m") = 50;
  parameter Real d(unit="mm") = 6;
  parameter Real Ptar(unit="bar")=320;
  parameter Modelica.SIunits.Force Fm = 0;
  Real x(unit="mm",start=0),v(unit="m/s",start=0);
  Modelica.SIunits.Force F;
  protected
  parameter Real x1(unit="mm")=0,x2(unit="mm")=3.3;
  parameter Real A(unit="mm2")=0.25*Modelica.Constants.pi*d^2;
  parameter Modelica.SIunits.Force F0 = 0.1*Ptar*A;
  parameter Real e = 0.1;
  Boolean impactL,canBounceL(start=true),impactR,canBounceR(start=true);
equation
  v = der(0.001*x);
  if canBounceL and canBounceR then
    m*der(v)=F;
  else
    m*der(v)=0;
  end if;
  F = Fm-c*v-k*x-F0;
algorithm
  impactL := x<=x1;
  impactR := x>=x2;
  when {impactL, impactL and v<=0} then
    if edge(impactL) then
      canBounceL := pre(v)<=0;
      reinit(v,-e*pre(v));
      reinit(x,x1);
    else
      reinit(v,0.0);
      canBounceL := false;
    end if;
  end when;
  when {impactR, impactR and v>=0} then
    if edge(impactR) then
      canBounceR := pre(v)>=0;
      reinit(v,-e*pre(v));
    else
      reinit(v,0.0);
      canBounceR := false;
    end if;
  end when;
end VS30CC;

The problem is, of course, with Fm = 1000 N works as before (x = 0.5 mm), even with Fm = 10000 works well (x = x2), but if imposed Fm = 0 the resulting chart is chilling: that the stroke instantly get to x = -8.5 mm before stopping at x = -1.4 mm.

Maybe I done something wrong. Do you know where is the mistake? Moreover, what is the best way to model efficiently rigid constraints? Is there any example of a body with two constraints that kind of my problem?

Thank you very much.
Best Regards,
Paolo Ferraresi, Italy.
fp.box@alice.it

There are 0 guests and 0 other users also viewing this topic
You are here: