475 lines
16 KiB
Python
475 lines
16 KiB
Python
from __future__ import annotations
|
|
|
|
import unittest
|
|
from math import exp
|
|
|
|
from app.simulation.components.amesim.flow.pipes import AmesimPnl0001
|
|
from app.simulation.components.amesim.media.mediums import (
|
|
AmesimHeliumPengRobinsonMedium,
|
|
)
|
|
from app.simulation.components.amesim.mechanical.pistons import AmesimPnrp17
|
|
from app.simulation.components.amesim.mechanical.translational import (
|
|
AmesimLstp00a,
|
|
AmesimMecmas21,
|
|
)
|
|
from app.simulation.components.amesim.storage.chambers import AmesimPnch012
|
|
from app.simulation.core.medium import IdealGasMedium
|
|
|
|
|
|
class AnalyticTangentPrimitiveTests(unittest.TestCase):
|
|
def assert_tangent_close(
|
|
self,
|
|
actual: float,
|
|
expected: float,
|
|
*,
|
|
relative: float = 5.0e-5,
|
|
absolute: float = 1.0e-8,
|
|
) -> None:
|
|
self.assertAlmostEqual(
|
|
actual,
|
|
expected,
|
|
delta=max(absolute, relative * max(abs(actual), abs(expected))),
|
|
)
|
|
|
|
def test_peng_robinson_m_u_v_linearization_matches_centered_difference(self) -> None:
|
|
medium = AmesimHeliumPengRobinsonMedium()
|
|
pressure = 15.3e6
|
|
temperature = 293.15
|
|
volume = 0.01
|
|
mass = medium.density(pressure, temperature) * volume
|
|
energy = mass * medium.specific_internal_energy_at_pressure(
|
|
pressure,
|
|
temperature,
|
|
)
|
|
linearization = medium.linearize_properties_from_mU(
|
|
mass,
|
|
energy,
|
|
volume,
|
|
(1.0, 0.0, 0.0),
|
|
(0.0, 1.0, 0.0),
|
|
(0.0, 0.0, 1.0),
|
|
)
|
|
self.assertTrue(linearization.valid, linearization.reason)
|
|
|
|
arguments = (mass, energy, volume)
|
|
steps = (mass * 1.0e-6, abs(energy) * 1.0e-6, volume * 1.0e-6)
|
|
for direction, step in enumerate(steps):
|
|
lower = list(arguments)
|
|
upper = list(arguments)
|
|
lower[direction] -= step
|
|
upper[direction] += step
|
|
lower_props = medium.properties_from_mU(*lower)
|
|
upper_props = medium.properties_from_mU(*upper)
|
|
for field in ("p", "T", "rho", "u", "h"):
|
|
finite_difference = (
|
|
getattr(upper_props, field) - getattr(lower_props, field)
|
|
) / (2.0 * step)
|
|
tangent = getattr(linearization.tangents, field)[direction]
|
|
self.assert_tangent_close(
|
|
tangent,
|
|
finite_difference,
|
|
relative=2.0e-4,
|
|
absolute=1.0e-6,
|
|
)
|
|
|
|
def test_pnrp17_geometry_and_force_tangent_is_exact(self) -> None:
|
|
piston = AmesimPnrp17("piston", IdealGasMedium(), dp=0.2, dr=0.01, x0=0.1)
|
|
piston.port_4.x = 0.02
|
|
piston.port_5.x = 0.05
|
|
piston.port_4.v = -0.2
|
|
piston.port_5.v = 0.3
|
|
piston.port_1.p = 2.0e5
|
|
result = piston.linearize_geometry_and_force(
|
|
(1.0, 0.0),
|
|
(0.0, 2.0),
|
|
(3.0, 0.0),
|
|
(0.0, 4.0),
|
|
(5.0, 6.0),
|
|
)
|
|
area = piston.effective_area
|
|
self.assertTrue(result.valid)
|
|
self.assertEqual(result.volume_tangent, (-area, 2.0 * area))
|
|
self.assertEqual(result.volume_flow_tangent, (-3.0 * area, 4.0 * area))
|
|
self.assertEqual(result.pressure_force_tangent, (5.0 * area, 6.0 * area))
|
|
|
|
def test_pnch012_balance_tangent_matches_directional_difference(self) -> None:
|
|
chamber = AmesimPnch012(
|
|
"chamber",
|
|
IdealGasMedium(),
|
|
cvol0=0.02,
|
|
kth=3.0,
|
|
sth=0.4,
|
|
p0=2.0e5,
|
|
T0=300.0,
|
|
)
|
|
chamber.port_1.volume = 0.003
|
|
chamber.port_1.volume_flow = 2.0e-4
|
|
flows = (0.02, -0.01, 0.005, -0.004)
|
|
enthalpies = {
|
|
"port_1": 330000.0,
|
|
"port_2": 310000.0,
|
|
"port_3": 320000.0,
|
|
"port_4": 300000.0,
|
|
}
|
|
for port_name, flow in zip(
|
|
("port_1", "port_2", "port_3", "port_4"),
|
|
flows,
|
|
strict=True,
|
|
):
|
|
chamber.get_port(port_name).m_flow = flow
|
|
|
|
dm = (1.0e-5,)
|
|
dU = (2.0,)
|
|
dV = (3.0e-6,)
|
|
dVdt = (-4.0e-5,)
|
|
dq = {
|
|
"port_1": (2.0e-3,),
|
|
"port_2": (-1.0e-3,),
|
|
"port_3": (3.0e-3,),
|
|
"port_4": (-2.0e-3,),
|
|
}
|
|
dh = {
|
|
"port_1": (20.0,),
|
|
"port_2": (-10.0,),
|
|
"port_3": (30.0,),
|
|
"port_4": (-20.0,),
|
|
}
|
|
result = chamber.linearize_state_derivative(
|
|
enthalpies,
|
|
state_mass_tangent=dm,
|
|
state_energy_tangent=dU,
|
|
external_volume_tangent=dV,
|
|
external_volume_rate_tangent=dVdt,
|
|
port_mass_flow_tangents=dq,
|
|
connected_h_tangents=dh,
|
|
)
|
|
self.assertTrue(result.valid, result.reason)
|
|
|
|
original_state = chamber.get_state_vector()
|
|
original_volume = chamber.port_1.volume
|
|
original_volume_flow = chamber.port_1.volume_flow
|
|
epsilon = 1.0e-4
|
|
|
|
def evaluate(sign: float) -> list[float]:
|
|
chamber.set_state_vector(
|
|
[
|
|
original_state[0] + sign * epsilon * dm[0],
|
|
original_state[1] + sign * epsilon * dU[0],
|
|
]
|
|
)
|
|
chamber.port_1.volume = original_volume + sign * epsilon * dV[0]
|
|
chamber.port_1.volume_flow = (
|
|
original_volume_flow + sign * epsilon * dVdt[0]
|
|
)
|
|
perturbed_h = {}
|
|
for port_name in enthalpies:
|
|
port = chamber.get_port(port_name)
|
|
base_flow = flows[int(port_name[-1]) - 1]
|
|
port.m_flow = base_flow + sign * epsilon * dq[port_name][0]
|
|
perturbed_h[port_name] = (
|
|
enthalpies[port_name] + sign * epsilon * dh[port_name][0]
|
|
)
|
|
return chamber.state_derivative_from_ports(perturbed_h)
|
|
|
|
lower = evaluate(-1.0)
|
|
upper = evaluate(1.0)
|
|
chamber.set_state_vector(original_state)
|
|
chamber.port_1.volume = original_volume
|
|
chamber.port_1.volume_flow = original_volume_flow
|
|
for port_name, flow in zip(
|
|
("port_1", "port_2", "port_3", "port_4"),
|
|
flows,
|
|
strict=True,
|
|
):
|
|
chamber.get_port(port_name).m_flow = flow
|
|
for row in range(2):
|
|
finite_difference = (upper[row] - lower[row]) / (2.0 * epsilon)
|
|
self.assert_tangent_close(
|
|
result.tangents[row][0],
|
|
finite_difference,
|
|
relative=1.0e-5,
|
|
absolute=1.0e-6,
|
|
)
|
|
|
|
def test_pnl0001_local_flow_slope_and_balance_tangent(self) -> None:
|
|
medium = IdealGasMedium()
|
|
pipe = AmesimPnl0001(
|
|
"pipe",
|
|
medium,
|
|
diam=0.02,
|
|
le=0.5,
|
|
kth=2.0,
|
|
p0=2.0e5,
|
|
T0=300.0,
|
|
)
|
|
slope = pipe.linearize_mass_flow(2.2e5, 1.8e5, 300.0)
|
|
self.assertTrue(slope.valid, slope.reason)
|
|
self.assertGreater(slope.partial_p_1, 0.0)
|
|
self.assertLess(slope.partial_p_2, 0.0)
|
|
self.assertLess(slope.partial_temperature, 0.0)
|
|
boundary = pipe.linearize_mass_flow(2.0e5, 2.0e5, 300.0)
|
|
self.assertFalse(boundary.valid)
|
|
self.assertEqual(boundary.reason, "flow_direction_boundary")
|
|
|
|
pipe.port_1.m_flow = 0.02
|
|
pipe.port_2.m_flow = -0.01
|
|
connected_h = {"port_1": 330000.0, "port_2": 310000.0}
|
|
dm = (1.0e-6,)
|
|
dU = (0.3,)
|
|
dq = {"port_1": (2.0e-3,), "port_2": (-1.0e-3,)}
|
|
dh = {"port_1": (20.0,), "port_2": (-10.0,)}
|
|
result = pipe.linearize_state_derivative(
|
|
connected_h,
|
|
state_mass_tangent=dm,
|
|
state_energy_tangent=dU,
|
|
port_mass_flow_tangents=dq,
|
|
connected_h_tangents=dh,
|
|
)
|
|
self.assertTrue(result.valid, result.reason)
|
|
|
|
original_state = pipe.get_state_vector()
|
|
epsilon = 1.0e-4
|
|
|
|
def evaluate(sign: float) -> list[float]:
|
|
pipe.set_state_vector(
|
|
[
|
|
original_state[0] + sign * epsilon * dm[0],
|
|
original_state[1] + sign * epsilon * dU[0],
|
|
]
|
|
)
|
|
pipe.port_1.m_flow = 0.02 + sign * epsilon * dq["port_1"][0]
|
|
pipe.port_2.m_flow = -0.01 + sign * epsilon * dq["port_2"][0]
|
|
return pipe.state_derivative_from_ports(
|
|
{
|
|
name: value + sign * epsilon * dh[name][0]
|
|
for name, value in connected_h.items()
|
|
}
|
|
)
|
|
|
|
lower = evaluate(-1.0)
|
|
upper = evaluate(1.0)
|
|
for row in range(2):
|
|
finite_difference = (upper[row] - lower[row]) / (2.0 * epsilon)
|
|
self.assert_tangent_close(
|
|
result.tangents[row][0],
|
|
finite_difference,
|
|
relative=1.0e-5,
|
|
absolute=1.0e-6,
|
|
)
|
|
|
|
def test_lstp_fixed_mode_tangent_and_boundary_guard(self) -> None:
|
|
contact = AmesimLstp00a(
|
|
"contact",
|
|
IdealGasMedium(),
|
|
gap0=0.001,
|
|
kcont=2000.0,
|
|
rcont=10.0,
|
|
Pdis=0.0005,
|
|
)
|
|
contact.port_1.x = 0.002
|
|
contact.port_2.x = 0.0
|
|
contact.port_1.v = 0.3
|
|
contact.port_2.v = -0.1
|
|
result = contact.linearize_contact_force(
|
|
(0.4,),
|
|
(-0.2,),
|
|
(0.3,),
|
|
(-0.1,),
|
|
)
|
|
self.assertTrue(result.valid, result.reason)
|
|
epsilon = 1.0e-7
|
|
originals = (
|
|
contact.port_1.x,
|
|
contact.port_2.x,
|
|
contact.port_1.v,
|
|
contact.port_2.v,
|
|
)
|
|
|
|
def evaluate(sign: float) -> float:
|
|
contact.port_1.x = originals[0] + sign * epsilon * 0.4
|
|
contact.port_2.x = originals[1] + sign * epsilon * -0.2
|
|
contact.port_1.v = originals[2] + sign * epsilon * 0.3
|
|
contact.port_2.v = originals[3] + sign * epsilon * -0.1
|
|
return contact.contact_force
|
|
|
|
finite_difference = (evaluate(1.0) - evaluate(-1.0)) / (2.0 * epsilon)
|
|
self.assert_tangent_close(result.force_tangent[0], finite_difference)
|
|
|
|
contact.clear_causal_contact()
|
|
contact.port_1.x = contact.gap0
|
|
contact.port_2.x = 0.0
|
|
boundary = contact.linearize_contact_force(
|
|
(1.0,),
|
|
(0.0,),
|
|
(0.0,),
|
|
(0.0,),
|
|
)
|
|
self.assertFalse(boundary.valid)
|
|
self.assertEqual(boundary.reason, "contact_mode_boundary")
|
|
|
|
def test_lstp_causal_contact_retains_a_differentiable_local_mode(self) -> None:
|
|
contact = AmesimLstp00a(
|
|
"causal_contact",
|
|
IdealGasMedium(),
|
|
gap0=0.0,
|
|
kcont=3000.0,
|
|
rcont=12.0,
|
|
Pdis=0.001,
|
|
)
|
|
contact.port_1.x = 1000.0
|
|
contact.port_2.x = 1000.0
|
|
contact.port_1.v = 0.2
|
|
contact.port_2.v = -0.1
|
|
penetration = 0.002
|
|
damping_fraction = 1.0 - exp(-penetration / contact.Pdis)
|
|
force = (
|
|
contact.kcont * penetration
|
|
+ damping_fraction * contact.rcont * contact.penetration_velocity
|
|
)
|
|
contact.set_causal_contact(penetration=penetration, force=force)
|
|
result = contact.linearize_contact_force(
|
|
(0.3,),
|
|
(-0.2,),
|
|
(0.4,),
|
|
(-0.1,),
|
|
)
|
|
self.assertTrue(result.valid, result.reason)
|
|
epsilon = 1.0e-6
|
|
originals = (
|
|
contact.port_1.x,
|
|
contact.port_2.x,
|
|
contact.port_1.v,
|
|
contact.port_2.v,
|
|
)
|
|
|
|
def evaluate(sign: float) -> float:
|
|
contact.port_1.x = originals[0] + sign * epsilon * 0.3
|
|
contact.port_2.x = originals[1] + sign * epsilon * -0.2
|
|
contact.port_1.v = originals[2] + sign * epsilon * 0.4
|
|
contact.port_2.v = originals[3] + sign * epsilon * -0.1
|
|
return contact.contact_force
|
|
|
|
finite_difference = (evaluate(1.0) - evaluate(-1.0)) / (2.0 * epsilon)
|
|
self.assert_tangent_close(
|
|
result.force_tangent[0],
|
|
finite_difference,
|
|
relative=2.0e-5,
|
|
)
|
|
|
|
def test_mecmas_soft_endstop_tangents_match_centered_difference(self) -> None:
|
|
cases = (
|
|
{
|
|
"name": "lower",
|
|
"x0": -0.02,
|
|
"xmin": 0.0,
|
|
"xmax": 1.0,
|
|
},
|
|
{
|
|
"name": "upper",
|
|
"x0": 1.02,
|
|
"xmin": 0.0,
|
|
"xmax": 1.0,
|
|
},
|
|
)
|
|
for case in cases:
|
|
with self.subTest(side=case["name"]):
|
|
mass = AmesimMecmas21(
|
|
f"mass_{case['name']}",
|
|
IdealGasMedium(),
|
|
mass=2.0,
|
|
useFriction=1.0,
|
|
stoptype=2.0,
|
|
discContactOption=1.0,
|
|
Kbmin=1000.0,
|
|
Dbmin=20.0,
|
|
Pdmin=0.1,
|
|
Kbmax=1200.0,
|
|
Dbmax=30.0,
|
|
Pdmax=0.1,
|
|
v0=0.3,
|
|
x0=case["x0"],
|
|
xmin=case["xmin"],
|
|
xmax=case["xmax"],
|
|
)
|
|
mass.port_1.f = 5.0
|
|
mass.port_2.f = -1.0
|
|
result = mass.linearize_state_derivative(
|
|
(0.2,),
|
|
(-0.1,),
|
|
(-0.25,),
|
|
(0.4,),
|
|
constraint_mode="free",
|
|
)
|
|
self.assertTrue(result.valid, result.reason)
|
|
epsilon = 1.0e-7
|
|
original_v = mass.v
|
|
original_x = mass.x
|
|
|
|
def evaluate(sign: float) -> float:
|
|
mass.v = original_v + sign * epsilon * -0.25
|
|
mass.x = original_x + sign * epsilon * 0.4
|
|
mass.port_1.f = 5.0 + sign * epsilon * 0.2
|
|
mass.port_2.f = -1.0 + sign * epsilon * -0.1
|
|
return mass.unconstrained_acceleration()
|
|
|
|
finite_difference = (
|
|
evaluate(1.0) - evaluate(-1.0)
|
|
) / (2.0 * epsilon)
|
|
self.assert_tangent_close(
|
|
result.tangents[0][0],
|
|
finite_difference,
|
|
relative=1.0e-6,
|
|
)
|
|
|
|
def test_mecmas_free_and_fixed_mode_tangents(self) -> None:
|
|
mass = AmesimMecmas21(
|
|
"mass",
|
|
IdealGasMedium(),
|
|
mass=2.0,
|
|
useFriction=2.0,
|
|
rvisc=3.0,
|
|
wind=0.5,
|
|
fcoul=0.0,
|
|
stoptype=4.0,
|
|
v0=2.0,
|
|
)
|
|
mass.port_1.f = 10.0
|
|
mass.port_2.f = -2.0
|
|
result = mass.linearize_state_derivative(
|
|
(0.3,),
|
|
(-0.1,),
|
|
(0.2,),
|
|
(0.4,),
|
|
constraint_mode="free",
|
|
)
|
|
self.assertTrue(result.valid, result.reason)
|
|
epsilon = 1.0e-6
|
|
original_v = mass.v
|
|
original_x = mass.x
|
|
|
|
def evaluate(sign: float) -> float:
|
|
mass.v = original_v + sign * epsilon * 0.2
|
|
mass.x = original_x + sign * epsilon * 0.4
|
|
mass.port_1.f = 10.0 + sign * epsilon * 0.3
|
|
mass.port_2.f = -2.0 + sign * epsilon * -0.1
|
|
return mass.unconstrained_acceleration()
|
|
|
|
finite_difference = (evaluate(1.0) - evaluate(-1.0)) / (2.0 * epsilon)
|
|
self.assert_tangent_close(result.tangents[0][0], finite_difference)
|
|
self.assertEqual(result.tangents[1], (0.2,))
|
|
|
|
mass.set_constraint_motion(0.0, velocity=0.0)
|
|
fixed = mass.linearize_state_derivative(
|
|
(1.0,),
|
|
(1.0,),
|
|
(1.0,),
|
|
(1.0,),
|
|
constraint_mode="lower",
|
|
)
|
|
self.assertTrue(fixed.valid, fixed.reason)
|
|
self.assertEqual(fixed.tangents, ((0.0,), (0.0,)))
|
|
|
|
|
|
if __name__ == "__main__":
|
|
unittest.main()
|