Added inertia mesure with pyoccad

twiinIT · Jul 31, 2023 · b168fb4 · b168fb4
1 parent 4789d80
commit b168fb4
Show file tree

Hide file tree

Showing 2 changed files with 23 additions and 0 deletions.
diff --git a/rocket_twin/systems/tank/tank.py b/rocket_twin/systems/tank/tank.py
@@ -1,4 +1,9 @@
+import numpy as np
 from cosapp.base import System
+from OCC.Core.BRepGProp import brepgprop
+from OCC.Core.gp import gp_Vec, gp_XYZ
+from OCC.Core.GProp import GProp_GProps
+from pyoccad.create import CreateCylinder
 
 
 class Tank(System):
@@ -24,9 +29,17 @@ class Tank(System):
     def setup(self):
 
         # Geometry
+        self.add_inward("radius", 1.0, desc="Tank base radius", unit="m")
+        self.add_inward("height", 1.0, desc="Tank height", unit="m")
         self.add_inward("weight_s", 1.0, desc="Structure weight", unit="kg")
         self.add_inward("weight_max", 5.0, desc="Maximum fuel capacity", unit="kg")
 
+        # pyoccad model
+        shape = CreateCylinder.from_base_and_dir(
+            np.zeros(3), gp_Vec(gp_XYZ(0, 0, self.height)), self.radius
+        )
+        self.add_inward("shape", shape, desc="pyoccad tank model")
+
         # Inputs
         self.add_inward("w_in", 0.0, desc="Fuel income rate", unit="kg/s")
 
@@ -38,10 +51,17 @@ def setup(self):
         self.add_outward("weight", 1.0, desc="Weight", unit="kg")
         self.add_outward("cg", 1.0, desc="Center of gravity", unit="m")
         self.add_outward("w_out", 0.0, desc="Fuel output rate", unit="kg/s")
+        self.add_outward("I", np.zeros((3, 3)), desc="Tank inertia tensor", unit="kg*m**2")
 
         # Transient
         self.add_transient("weight_p", der="w_in - w_out", desc="Propellant weight")
 
     def compute(self):
         self.w_out = self.w_out_max * self.w_command
         self.weight = self.weight_s + self.weight_p
+
+        v_prop = GProp_GProps()
+        brepgprop.VolumeProperties(self.shape, v_prop)
+        inertia = v_prop.MatrixOfInertia()
+        for i, j in zip(range(3), range(3)):
+            self.I[i, j] = inertia.Value(i + 1, j + 1)
diff --git a/rocket_twin/tests/test_tank.py b/rocket_twin/tests/test_tank.py
@@ -19,6 +19,9 @@ def test_fuel(self):
         sys.run_drivers()
 
         np.testing.assert_allclose(sys.weight, 16.0, atol=10 ** (-10))
+        np.testing.assert_allclose(
+            sys.I.diagonal(), np.array([np.pi / 3, np.pi / 3, np.pi / 2]), atol=10 ** (-10)
+        )
 
     def test_flight(self):
         sys = Tank("sys")