I’m trying to model the 2 Body Problem in OpenModelica (as a first step to model the 3 Body Problem) but I’m having some problems.
I’m using the LineForceWithMass object in order to have the force between the 2 objects on the line connecting the 2 but when i do it the position of the 2 bodies isn’t update even if the force is computed in the correct way.

I realized a class “AttractionForce” that computes the force between the 2 bodies:

model AttractionForce
  
  import Modelica.Mechanics.Translational.Interfaces;
  extends Modelica.Mechanics.Translational.Interfaces.PartialCompliant;
  
  constant Real G = 6.67430e-11 "Gravitational constant";
  parameter SI.Mass m1;
  parameter SI.Mass m2;
  
  equation
  f = (G*m1*m2)/(s_rel^2);

annotation(
    Icon(graphics = {Polygon(origin = {10, 0}, lineColor = {255, 85, 0}, fillColor = {255, 85, 0}, fillPattern = FillPattern.Solid, points = {{-70, 20}, {-70, -20}, {30, -20}, {30, -40}, {70, 0}, {30, 40}, {30, 20}, {-70, 20}, {-70, 20}}), Polygon(origin = {-60, 0}, lineColor = {255, 85, 0}, fillColor = {255, 85, 0}, fillPattern = FillPattern.Solid, points = {{-20, 0}, {20, 40}, {20, -40}, {-20, 0}, {-20, 0}})}));

end AttractionForce;

Then I used this class with the LineForceWithMass object to model the gravitational attraction between the 2 bodies:

model GravitationalForce
  
  import Modelica.Mechanics.MultiBody.Interfaces;
  import Modelica.Mechanics.MultiBody.Forces;
  import Modelica.Mechanics.MultiBody.Frames;
  extends Interfaces.PartialTwoFrames;
  
  parameter Boolean animation = true;
  parameter Boolean showMass = true;
  parameter SI.Distance s_small = 1e-10;
  parameter SI.Mass m1;
  parameter SI.Mass m2;
 
  Modelica.Mechanics.MultiBody.Forces.LineForceWithMass lineForce(animateLine = animation, animateMass = showMass, fixedRotationAtFrame_a = true, fixedRotationAtFrame_b = true, m = 0, s_small = s_small) annotation(
    Placement(transformation(extent = {{-20, -20}, {20, 20}})));
 
  SI.Position r_rel_a[3] "Position vector from origin of frame_a to origin of frame_b, resolved in frame_a";
  Real e_a[3](each final unit = "1") "Unit vector on the line connecting the origin of frame_a with the origin of frame_b resolved in frame_a (directed from frame_a to frame_b)";
  SI.Force f "Line force acting on frame_a and on frame_b (positive, if acting on frame_b and directed from frame_a to frame_b)";
  SI.Distance length "Distance between the origin of frame_a and the origin of frame_b";
  SI.Position s "(Guarded) distance between the origin of frame_a and the origin of frame_b (>= s_small))";
  SI.Position r_rel_0[3] "Position vector from frame_a to frame_b resolved in world frame";
  Real e_rel_0[3](each final unit = "1") "Unit vector in direction from frame_a to frame_b, resolved in world frame";
 
  AttractionForce attractionForce( m1 = m1, m2 = m2)  annotation(
    Placement(visible = true, transformation(origin = {0, 44}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));

equation
  r_rel_a = Frames.resolve2(frame_a.R, r_rel_0);
  e_a = r_rel_a/s;
  f = attractionForce.f;
  length = lineForce.length;
  s = lineForce.s;
  r_rel_0 = lineForce.frame_a.r_0 - lineForce.frame_b.r_0;
  e_rel_0 = r_rel_0/s;
 
  connect(frame_a, lineForce.frame_a) annotation(
    Line(points = {{-100, 0}, {-20, 0}}));
  connect(lineForce.frame_b, frame_b) annotation(
    Line(points = {{20, 0}, {100, 0}}, color = {95, 95, 95}));
  connect(lineForce.flange_a, attractionForce.flange_a) annotation(
    Line(points = {{-12, 20}, {-10, 20}, {-10, 44}}, color = {0, 127, 0}));
  connect(lineForce.flange_b, attractionForce.flange_b) annotation(
    Line(points = {{12, 20}, {10, 20}, {10, 44}}, color = {0, 127, 0}));


annotation(
    Icon(graphics = {Polygon(origin = {-50, 0}, lineColor = {0, 0, 255}, fillColor = {0, 0, 255}, fillPattern = FillPattern.Solid, points = {{-30, 0}, {-10, 40}, {-10, 20}, {28, 20}, {30, 20}, {30, -20}, {-10, -20}, {-10, -40}, {-30, 0}, {-30, 0}}), Polygon(origin = {50, 0}, lineColor = {0, 0, 255}, fillColor = {0, 0, 255}, fillPattern = FillPattern.Solid, points = {{-30, 20}, {-30, -18}, {-30, -20}, {10, -20}, {10, -40}, {30, 0}, {10, 40}, {10, 20}, {8, 20}, {-30, 20}})}));
end GravitationalForce;

Then i just put everything togheter to simulate the problem but it doesn’t work and i can’t manage to understand why (it’s one of my first times using OpenModelica so I don’t know if there’s something obvious that I’m overlooking).

model TwoBodyProblem
  
  import Modelica.Math.Vectors.norm;
  import Modelica.SIunits;
  
  parameter SI.Mass m1 = 1e24 "Mass of first body";
  parameter SI.Mass m2 = 1e24 "Mass of second body";
  
  inner Modelica.Mechanics.MultiBody.World world(animateGravity = true, axisDiameter = 1e8, axisLength = 1e9, gravitySphereDiameter = 0.1, gravityType = 
   Modelica.Mechanics.MultiBody.Types.GravityTypes.NoGravity, mu = 0) annotation(
    Placement(visible = true, transformation(origin = {16, -54}, extent = {{-80, -20}, {-60, 0}}, rotation = 0)));
 
  Modelica.Mechanics.MultiBody.Parts.Body planet1(I_11 = 0.1, I_22 = 0.1, I_33 = 0.1, m = m1, r_0(start = {1e10, 0, 0}), r_CM = {0, 0, 0}, sphereDiameter = 1e9)  annotation(
    Placement(visible = true, transformation(origin = {-60, 40}, extent = {{-10, -10}, {10, 10}}, rotation = 180)));
  Modelica.Mechanics.MultiBody.Parts.Body planet2(I_11 = 0.1, I_22 = 0.1, I_33 = 0.1, m = m2, r_0(start = {0, 1e10, 0}), r_CM = {0, 0, 0}, sphereDiameter = 1e9)  annotation(
    Placement(visible = true, transformation(origin = {36, 40}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
 GravitationalForce gravitationalForce(m1 = m1, m2 = m2) annotation(
    Placement(visible = true, transformation(origin = {-14, 40}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
equation
  connect(planet2.frame_a, gravitationalForce.frame_b) annotation(
    Line(points = {{26, 40}, {-4, 40}}, color = {95, 95, 95}));
 connect(planet1.frame_a, gravitationalForce.frame_a) annotation(
    Line(points = {{-50, 40}, {-24, 40}}, color = {95, 95, 95}));
  annotation(
    Icon(graphics = {Ellipse(origin = {-43, 38}, lineColor = {170, 85, 255}, fillColor = {170, 85, 255}, fillPattern = FillPattern.Solid, extent = {{-27, 28}, {27, -28}}), Ellipse(origin = {54, -55}, lineColor = {85, 255, 0}, fillColor = {0, 255, 0}, fillPattern = FillPattern.Solid, extent = {{-34, 35}, {34, -35}}), Ellipse(origin = {7, -12}, extent = {{-71, 68}, {71, -68}})}));
end TwoBodyProblem;