i've run into a problem with out of domain error involving modelica trigonometric functions and i just dont know what i could do.
My model is that of a nonlinear rotational spring (torque dependant on geometric relations) + a hard stop (using standard backlash element).
Model schematic with elastogap - doesnt work
The error i get in OpenModelica is:
acos(variable) outside domain -1.0 <= -1 <= +1.0
That means the error appears to be in my custom spring model.
Now, the ordinary solution i thought should work is to limit the variable so acos stays in domain using if-conditions. I also followed the advice given here about the noEvent() operator. This did not work however. I even went overboard with the "safekeeping", to no avail.
For completions sake, this is the geometric relation of the spring. I'm positive that it is correct.
I use an input ramp from phi= -2° to 108°. The elastogap allows force-free movement from phi = 0° to 106°.
Without the backlash element my spring works as i expect . No issues with out of domain. Together with this backlash i get the domain issue, even though i added "domain-limits" before every trigonometric function in my custom spring model.
Does anyone have an idea what the issue could be/ what i could try?
I noticed that if i change the range of the backlash to have different limits (for example b=105.9*pi/180 ; phi_rel0=-52.95*pi/180) it also works without domain error, but i would like to get to the bottom of this issue.
The model (with the error):
package test
model LD_Abwickel_ZF_VAL
Modelica.Blocks.Sources.Ramp ramp2(duration = 1.24, height = 122, offset = -2);
Modelica.Mechanics.Rotational.Components.ElastoBacklash2 elastoBacklash(b = 1.85005, c = 1e4, d = 1e2, phi_rel0 = -0.925025);
Modelica.Mechanics.Rotational.Components.Fixed fixed2;
Modelica.Mechanics.Rotational.Sources.Position position1(exact = true);
Modelica.Blocks.Math.UnitConversions.From_deg from_deg1;
test.LD_Pendelarm_ZugfederF2 lD_Pendelarm_ZugfederF21;
equation
connect(position1.flange, lD_Pendelarm_ZugfederF21.flange);
connect(elastoBacklash.flange_a, lD_Pendelarm_ZugfederF21.flange);
connect(ramp2.y, from_deg1.u);
connect(from_deg1.y, position1.phi_ref);
connect(elastoBacklash.flange_b, fixed2.flange);
end LD_Abwickel_ZF_VAL;
model LD_Pendelarm_ZugfederF2
Modelica.Mechanics.Rotational.Interfaces.Flange_a flange;
.Modelica.SIunits.Angle phi(displayUnit = "deg");
.Modelica.SIunits.Torque M;
parameter .Modelica.SIunits.Length a(displayUnit = "mm") = 25.49e-3;
parameter .Modelica.SIunits.Length b(displayUnit = "mm") = 23.38e-3;
parameter .Modelica.SIunits.Length d(displayUnit = "mm") = 43.89e-3;
parameter .Modelica.SIunits.Length L0(displayUnit = "mm") = 59.5e-3 "Ungespannte Laenge d. Feder";
parameter .Modelica.SIunits.Length Lk(displayUnit = "mm") = 47.19e-3 "Laenge unbelasteter Federkörper";
parameter .Modelica.SIunits.Force F0 = 4.1 "innere Vorspannkraft Feder";
parameter .Modelica.SIunits.TranslationalSpringConstant R(displayUnit = "N/mm") = 0.868e+3 "Federkonstante";
parameter Integer z = 2 "Anzahl Federn";
.Modelica.SIunits.Length c(displayUnit = "mm");
.Modelica.SIunits.Angle delta(displayUnit = "deg", min = -2 * .Modelica.Constants.pi, max = 2 * .Modelica.Constants.pi, nominal = 0);
.Modelica.SIunits.Angle alpha(displayUnit = "deg", min = -2 * .Modelica.Constants.pi, max = 2 * .Modelica.Constants.pi, nominal = 0.5 * .Modelica.Constants.pi);
.Modelica.SIunits.Angle beta(displayUnit = "deg", min = -4 * .Modelica.Constants.pi, max = 4 * .Modelica.Constants.pi, nominal = .Modelica.Constants.pi);
.Modelica.SIunits.Angle psi(displayUnit = "deg", min = -4 * .Modelica.Constants.pi, max = 4 * .Modelica.Constants.pi, nominal = 0);
.Modelica.SIunits.Length L_gerade(displayUnit = "mm");
.Modelica.SIunits.Length L_bogen(displayUnit = "mm") = 46.1e-3;
.Modelica.SIunits.Length L_c(displayUnit = "mm") "rel. Federlaenge (gesamt)";
.Modelica.SIunits.Force F_c "Summe Federkraft";
Real cosalpha(min = -1, max = 1, nominal = 0);
Real cosdelta(min = -1, max = 1, nominal = 0);
Real cosbeta(min = -1, max = 1, nominal = 0);
equation
flange.phi = phi;
flange.tau = M;
if noEvent(phi <= 0) then
c = a - b;
cosdelta = 1;
delta = 0;
L_gerade = sqrt(c ^ 2 + d ^ 2);
cosalpha = -1;
alpha = .Modelica.Constants.pi;
else
c = sqrt(a ^ 2 + b ^ 2 - 2 * a * b * cos(phi));
cosdelta = (b ^ 2 - a ^ 2 - c ^ 2) / (-2 * a * c);
if noEvent(cosdelta <= 0.999999 and cosdelta >= (-0.999999)) then
cos(delta) = cosdelta;
else
delta = if noEvent(sign(cosdelta) > 0) then 0 else .Modelica.Constants.pi;
end if;
L_gerade = sqrt(c ^ 2 + d ^ 2 - 2 * c * d * cos(delta + .Modelica.Constants.pi * 0.5));
cosalpha = (a ^ 2 - b ^ 2 - c ^ 2) / (-2 * b * c);
if noEvent(cosalpha <= 0.999999 and cosalpha >= (-0.999999)) then
cos(alpha) = cosalpha;
else
alpha = if noEvent(sign(cosalpha) > 0) then 0 else .Modelica.Constants.pi;
end if;
end if;
cosbeta = (d ^ 2 - c ^ 2 - L_gerade ^ 2) / (-2 * c * L_gerade);
if noEvent(cosbeta <= 0.999999 and cosbeta >= (-0.999999)) then
cos(beta) = cosbeta;
else
beta = if noEvent(sign(cosbeta) > 0) then 0 else .Modelica.Constants.pi;
end if;
psi = alpha - beta;
L_c = L_gerade + L_bogen - (L0 - Lk) - Lk;
F_c = z * (R * L_c + F0);
M = F_c * sin(psi);
end LD_Pendelarm_ZugfederF2;
end test;
So of course it works fine on my machine, but from the code I know it is problematic and might fail randomly. The problem is that you have things like cos(beta) = cosbeta which is solved as beta = acos(cosbeta).
But OpenModelica added a bad equation for common subexpression elimination someTmpVar = acos(cosbeta) which always runs regardless of the if-equations you use to guard this. The expansion seems to happen Feels like a bug to report to OpenModelica since it happens during code generation and just produces lD_Pendelarm_ZugfederF21._delta = if noEvent(from_deg1.y <= 0.0) then 0.0 else if noEvent(lD_Pendelarm_ZugfederF21.cosdelta <= 0.999999) and noEvent(lD_Pendelarm_ZugfederF21.cosdelta >= -0.999999) then if noEvent(abs($TMP$VAR$70$0x1$COS - $TMP$VAR$70$0PREX$COS) < abs($TMP$VAR$70$0x2$COS - $TMP$VAR$70$0PREX$COS)) then $TMP$VAR$70$0x1$COS else $TMP$VAR$70$0x2$COS else if noEvent(sign(lD_Pendelarm_ZugfederF21.cosdelta) > 0) then 0.0 else 3.141592653589793, which won't run faster than before and might produce domain errors.
A probable workaround would be to use an algorithm section as it is not manipulated as heavily by OpenModelica.