diff --git a/PythonModels/components/amesim_mechanical.py b/PythonModels/components/amesim_mechanical.py index 1669d22..2f3c4a4 100644 --- a/PythonModels/components/amesim_mechanical.py +++ b/PythonModels/components/amesim_mechanical.py @@ -125,3 +125,44 @@ class AmesimMassFrictionEndstops: def viscous_friction_force(self, velocity_m_s: float) -> float: return -self.viscous_friction_n_per_m_per_s * velocity_m_s + def windage_force(self, velocity_m_s: float) -> float: + return -self.windage_n_per_m2_per_s2 * velocity_m_s * abs(velocity_m_s) + + def dry_friction_force(self, velocity_m_s: float) -> float: + if velocity_m_s > 0.0: + return -self.coulomb_friction_n + if velocity_m_s < 0.0: + return self.coulomb_friction_n + return 0.0 + + def limit_contact_force(self, displacement_m: float, velocity_m_s: float) -> float: + lower_force = self.lower_static_force_magnitude(displacement_m) + if lower_force > 0.0: + lower_force += max(-self.lower_damping_n_per_m_per_s * velocity_m_s, 0.0) + + upper_force = self.upper_static_force_magnitude(displacement_m) + if upper_force > 0.0: + upper_force += max(self.upper_damping_n_per_m_per_s * velocity_m_s, 0.0) + + return lower_force - upper_force + + def derivatives( + self, + *, + velocity_m_s: float, + displacement_m: float, + port_1_force_n: float = 0.0, + port_2_force_n: float = 0.0, + external_force_n: float = 0.0, + ) -> tuple[float, float]: + total_force = ( + port_1_force_n + + port_2_force_n + + external_force_n + + self.viscous_friction_force(velocity_m_s) + + self.windage_force(velocity_m_s) + + self.dry_friction_force(velocity_m_s) + + self.limit_contact_force(displacement_m, velocity_m_s) + ) + return total_force / self.mass_kg, velocity_m_s + diff --git a/PythonModels/systems/test_mql.py b/PythonModels/systems/test_mql.py index 3aad4a2..e91528f 100644 --- a/PythonModels/systems/test_mql.py +++ b/PythonModels/systems/test_mql.py @@ -3865,6 +3865,7 @@ class TestMqlSystem: self.pnl00r_assembly = self._build_pnl00r_assembly() self.node3_assembly = self._build_node3_assembly() self.node4_assembly = self._build_node4_assembly() + self.mechanical_assembly = self._build_mechanical_assembly() self.pneumatic_assembly = self._build_pneumatic_assembly() pneumatic_components = self._pneumatic_components_by_alias() for spec in COMPONENT_SPECS: @@ -3887,6 +3888,13 @@ class TestMqlSystem: return build_test_mql_pneumatic_assembly() + @staticmethod + def _build_mechanical_assembly(): + # Local import avoids a module cycle: test_mql_config reads constants from this module. + from PythonModels.systems.test_mql_mechanical import build_test_mql_mechanical_assembly + + return build_test_mql_mechanical_assembly() + def _build_pnl0001_assembly(self): from PythonModels.systems.test_mql_pneumatic_lines import ( build_test_mql_pnl0001_assembly, @@ -4112,6 +4120,27 @@ class TestMqlSystem: def pneumatic_state_vector(self) -> list[float]: return self.network.initial_state_vector() + def mechanical_mass_closure(self): + from PythonModels.systems.test_mql_mechanical import TestMqlMechanicalMassClosure + + return TestMqlMechanicalMassClosure(self.mechanical_assembly) + + def mechanical_state_vector(self) -> list[float]: + return self.mechanical_mass_closure().initial_state_vector() + + def simulate_mechanical_masses( + self, + config: SolveIVPConfig | None = None, + t_eval: list[float] | None = None, + ): + closure = self.mechanical_mass_closure() + return integrate_ode( + rhs=lambda t, state: closure.rhs(state), + initial_state=closure.initial_state_vector(), + config=config or SolveIVPConfig(), + t_eval=t_eval, + ) + @staticmethod def pneumatic_branch_spec( *, diff --git a/PythonModels/systems/test_mql_mechanical.py b/PythonModels/systems/test_mql_mechanical.py index 9c96fca..a8f43ce 100644 --- a/PythonModels/systems/test_mql_mechanical.py +++ b/PythonModels/systems/test_mql_mechanical.py @@ -63,6 +63,8 @@ class TestMqlMassEndstopSpec: stribeck_constant_m_s: float use_friction: bool stop_type: int + initial_velocity_m_s: float + initial_displacement_m: float data_paths: tuple[str, ...] def endstop(self) -> AmesimMassFrictionEndstops: @@ -141,6 +143,64 @@ class TestMqlMechanicalAssembly: ) +@dataclass(frozen=True) +class TestMqlMechanicalMassState: + alias: str + velocity_m_s: float + displacement_m: float + + def as_vector(self) -> list[float]: + return [self.velocity_m_s, self.displacement_m] + + +@dataclass(frozen=True) +class TestMqlMechanicalMassSnapshot: + states: tuple[TestMqlMechanicalMassState, ...] + + @property + def state_count(self) -> int: + return 2 * len(self.states) + + +class TestMqlMechanicalMassClosure: + def __init__(self, assembly: TestMqlMechanicalAssembly) -> None: + self.assembly = assembly + self.mass_aliases = tuple(assembly.masses) + + def initial_state_vector(self) -> list[float]: + state: list[float] = [] + for alias in self.mass_aliases: + spec = self.assembly.masses[alias] + state.extend([spec.initial_velocity_m_s, spec.initial_displacement_m]) + return state + + def snapshot(self, state_vector: list[float] | None = None) -> TestMqlMechanicalMassSnapshot: + values = self.initial_state_vector() if state_vector is None else list(state_vector) + if len(values) != 2 * len(self.mass_aliases): + raise ValueError("mechanical mass state vector requires two values per mass") + states = tuple( + TestMqlMechanicalMassState( + alias=alias, + velocity_m_s=values[2 * index], + displacement_m=values[2 * index + 1], + ) + for index, alias in enumerate(self.mass_aliases) + ) + return TestMqlMechanicalMassSnapshot(states=states) + + def rhs(self, state_vector: list[float]) -> list[float]: + snapshot = self.snapshot(state_vector) + derivatives: list[float] = [] + for state in snapshot.states: + mass = self.assembly.masses[state.alias].endstop() + acceleration, velocity = mass.derivatives( + velocity_m_s=state.velocity_m_s, + displacement_m=state.displacement_m, + ) + derivatives.extend([acceleration, velocity]) + return derivatives + + def build_test_mql_mechanical_assembly( config: TestMqlConfig | None = None, amesim_results: AmesimResults | None = None, @@ -155,7 +215,7 @@ def build_test_mql_mechanical_assembly( for component in config.components_by_submodel("PNRP17") } masses = { - component.alias: _build_mass(component, variable_catalog) + component.alias: _build_mass(component, variable_catalog, amesim_results) for component in config.components_by_submodel("MECMAS21") } elastic_endstops = { @@ -202,6 +262,7 @@ def _build_piston( def _build_mass( component: TestMqlResolvedComponent, variable_catalog: TestMqlVariableCatalog | None, + amesim_results: AmesimResults | None, ) -> TestMqlMassEndstopSpec: return TestMqlMassEndstopSpec( alias=component.alias, @@ -224,6 +285,8 @@ def _build_mass( stribeck_constant_m_s=component.parameter_value("astrib"), use_friction=bool(int(component.parameter_value("useFriction"))), stop_type=int(component.parameter_value("stoptype")), + initial_velocity_m_s=_initial_value(amesim_results, f"v1@{component.alias}"), + initial_displacement_m=_initial_value(amesim_results, f"x1@{component.alias}"), data_paths=_data_paths(variable_catalog, component.alias), ) @@ -263,6 +326,12 @@ def n_per_mm_per_s_to_n_per_m_per_s(value: float) -> float: return value * N_PER_MM_PER_S_TO_N_PER_M_PER_S +def _initial_value(amesim_results: AmesimResults | None, data_path: str) -> float: + if amesim_results is None: + return 0.0 + return float(amesim_results.series(data_path)[0]) + + def _data_paths( variable_catalog: TestMqlVariableCatalog | None, alias: str, diff --git a/tests/test_test_mql_mechanical.py b/tests/test_test_mql_mechanical.py index e0268bf..9831a73 100644 --- a/tests/test_test_mql_mechanical.py +++ b/tests/test_test_mql_mechanical.py @@ -4,8 +4,11 @@ import unittest from pathlib import Path from PythonModels.reporting.amesim_results import load_test_mql_amesim_results +from PythonModels.core.solver import SolveIVPConfig from PythonModels.reporting.test_mql_variables import build_test_mql_variable_catalog +from PythonModels.systems.test_mql import TestMqlSystem from PythonModels.systems.test_mql_mechanical import ( + TestMqlMechanicalMassClosure, build_test_mql_mechanical_assembly, circular_area, mm_to_m, @@ -37,6 +40,43 @@ class TestMqlMechanicalAssemblyTests(unittest.TestCase): self.assertEqual(len(self.assembly.force_connectors), 2) self.assertEqual(self.assembly.component_count, 46) + def test_mechanical_mass_closure_exposes_all_mecmas21_states(self) -> None: + closure = TestMqlMechanicalMassClosure(self.assembly) + state = closure.initial_state_vector() + snapshot = closure.snapshot(state) + rhs = closure.rhs(state) + + self.assertEqual(closure.mass_aliases, tuple(self.assembly.masses)) + self.assertEqual(len(closure.mass_aliases), 10) + self.assertEqual(len(state), 20) + self.assertEqual(snapshot.state_count, 20) + self.assertEqual(len(rhs), 20) + self.assertEqual(snapshot.states[0].alias, "mass_friction_endstops_10") + self.assertAlmostEqual( + snapshot.states[0].velocity_m_s, + self.amesim_results.series("v1@mass_friction_endstops_10")[0], + ) + self.assertAlmostEqual( + snapshot.states[0].displacement_m, + self.amesim_results.series("x1@mass_friction_endstops_10")[0], + ) + for index, mass_state in enumerate(snapshot.states): + self.assertAlmostEqual(rhs[2 * index + 1], mass_state.velocity_m_s) + + def test_system_simulates_mechanical_mass_state_closure(self) -> None: + system = TestMqlSystem() + + solution = system.simulate_mechanical_masses( + config=SolveIVPConfig(t_start=0.0, t_stop=1.0e-5, max_step=1.0e-6), + t_eval=[0.0, 1.0e-5], + ) + + self.assertTrue(solution.success) + self.assertEqual(len(system.mechanical_state_vector()), 20) + self.assertEqual(len(solution.t), 2) + self.assertEqual(len(solution.y), 20) + self.assertEqual([len(row) for row in solution.y], [2] * 20) + def test_piston_geometry_is_converted_to_si_units(self) -> None: piston = self.assembly.pistons["pn_brp2_8"]