"""Driven two-bar slider-crank flexible multibody validation journal.""" import sys import numpy as np sys.path.insert(0, "./server/") from mesh_generators.structured import build_structured_mesh from pytorsocae import TorsoCAESession CRANK_LENGTH = 0.25 ROD_LENGTH = 0.65 LINK_WIDTH = 0.05 LINK_THICKNESS = 0.04 PISTON_LENGTH = 0.12 PISTON_HEIGHT = 0.10 PISTON_THICKNESS = 0.08 INITIAL_CRANK_ANGLE = np.radians(30.0) YOUNG_MODULUS = 210.0e9 POISSON_RATIO = 0.3 DENSITY = 7800.0 DRIVE_TORQUE = 12.0 DRIVE_TRACTION = ( DRIVE_TORQUE / (CRANK_LENGTH * LINK_WIDTH * 0.5 * CRANK_LENGTH) ) TIME_STEP = 0.005 NUM_STEPS = 185 def _oriented_link(length, angle, origin, surface_base, volume_tag, divisions): mesh = build_structured_mesh( bounds=( (0.0, length), (-LINK_WIDTH / 2.0, LINK_WIDTH / 2.0), (-LINK_THICKNESS / 2.0, LINK_THICKNESS / 2.0), ), divisions=(divisions, 1, 1), surface_tags={ "xmin": surface_base + 1, "xmax": surface_base + 2, "ymin": surface_base + 3, "ymax": surface_base + 4, "zmin": surface_base + 5, "zmax": surface_base + 6, }, volume_tag=volume_tag, ) cosine, sine = np.cos(angle), np.sin(angle) rotation = np.array( ( (cosine, -sine, 0.0), (sine, cosine, 0.0), (0.0, 0.0, 1.0), ) ) mesh["x"] = np.asarray(mesh["x"], dtype=np.float64) @ rotation.T + origin return mesh def _combine_disconnected_meshes(*meshes): coordinates = [] topology = {} node_offset = 0 type_props = {} for mesh in meshes: points = np.asarray(mesh["x"], dtype=np.float64) coordinates.append(points) for element_type, block in mesh["topologies"].items(): destination = topology.setdefault( int(element_type), {"topology": [], "cell_data": []}, ) destination["topology"].append( np.asarray(block["topology"], dtype=np.int64) + node_offset ) destination["cell_data"].append( np.asarray(block["cell_data"], dtype=np.int32) ) type_props.update( {int(key): value for key, value in mesh["type_props"].items()} ) node_offset += len(points) for block in topology.values(): block["topology"] = np.vstack(block["topology"]) block["cell_data"] = np.concatenate(block["cell_data"]) return { "x": np.vstack(coordinates), "topologies": topology, "type_props": type_props, "gdim": 3, } def mechanism_geometry(): crank_pin = np.array( ( CRANK_LENGTH * np.cos(INITIAL_CRANK_ANGLE), CRANK_LENGTH * np.sin(INITIAL_CRANK_ANGLE), 0.0, ) ) piston_x = ( crank_pin[0] + np.sqrt(ROD_LENGTH**2 - crank_pin[1] ** 2) ) piston_pin = np.array((piston_x, 0.0, 0.0)) rod_angle = np.arctan2( piston_pin[1] - crank_pin[1], piston_pin[0] - crank_pin[0], ) return crank_pin, piston_pin, rod_angle def build_session() -> TorsoCAESession: """Build the driven slider-crank without executing it.""" crank_pin, piston_pin, rod_angle = mechanism_geometry() crank = _oriented_link( CRANK_LENGTH, INITIAL_CRANK_ANGLE, np.zeros(3), 100, 1, 3, ) connecting_rod = _oriented_link( ROD_LENGTH, rod_angle, crank_pin, 200, 2, 5, ) piston = build_structured_mesh( bounds=( ( piston_pin[0] - PISTON_LENGTH / 2.0, piston_pin[0] + PISTON_LENGTH / 2.0, ), (-PISTON_HEIGHT / 2.0, PISTON_HEIGHT / 2.0), (-PISTON_THICKNESS / 2.0, PISTON_THICKNESS / 2.0), ), divisions=(1, 1, 1), surface_tags={ "xmin": 301, "xmax": 302, "ymin": 303, "ymax": 304, "zmin": 305, "zmax": 306, }, volume_tag=3, ) mesh = _combine_disconnected_meshes(crank, connecting_rod, piston) session = TorsoCAESession() session.inline_mesh( mesh, mesh_id="slider_crank_mesh", name="Flexible Slider-Crank", volume_names={ 1: "Crank", 2: "Connecting Rod", 3: "Piston", }, ) for solid_name in ("Crank", "Connecting Rod", "Piston"): session.solid(solid_name).material( E=YOUNG_MODULUS, nu=POISSON_RATIO, rho=DENSITY, ) session.set_flexible_multibody( frames=[ { "id": "ground", "origin": [0.0, 0.0, 0.0], "grounded": True, }, {"id": "crank_root", "origin": [0.0, 0.0, 0.0]}, {"id": "crank_pin", "origin": crank_pin.tolist()}, {"id": "rod_root", "origin": crank_pin.tolist()}, {"id": "rod_tip", "origin": piston_pin.tolist()}, {"id": "piston", "origin": piston_pin.tolist()}, ], bodies=[ { "id": "crank_body", "kind": "flexible", "reference_frame_id": "crank_root", }, { "id": "rod_body", "kind": "flexible", "reference_frame_id": "rod_root", }, { "id": "piston_body", "kind": "flexible", "reference_frame_id": "piston", }, ], attachments=[ { "frame_id": "crank_root", "surface_tag": 101, "body_id": "crank_body", }, { "frame_id": "crank_pin", "surface_tag": 102, "body_id": "crank_body", }, { "frame_id": "rod_root", "surface_tag": 201, "body_id": "rod_body", }, { "frame_id": "rod_tip", "surface_tag": 202, "body_id": "rod_body", }, { "frame_id": "piston", "surface_tag": 301, "body_id": "piston_body", }, ], joints=[ { "id": "main_bearing", "frame_a": "ground", "frame_b": "crank_root", "kind": "revolute", "axis": "z", }, { "id": "crank_pin", "frame_a": "crank_pin", "frame_b": "rod_root", "kind": "revolute", "axis": "z", }, { "id": "wrist_pin", "frame_a": "rod_tip", "frame_b": "piston", "kind": "revolute", "axis": "z", }, { "id": "piston_guide", "frame_a": "ground", "frame_b": "piston", "free_dofs": ["tx", "tz", "rx", "ry"], }, ], kkt_augmentation=0.0, projection={ "tolerance": 1.0e-9, "max_iterations": 20, }, ) drive_direction = np.array( (-np.sin(INITIAL_CRANK_ANGLE), np.cos(INITIAL_CRANK_ANGLE), 0.0) ) session.surface(106).bc( "traction", values=(DRIVE_TRACTION * drive_direction).tolist(), follower=True, ) session.set_physics( "structural", submodel="flexible_multibody_dynamics", backend="dolfinx", ) session.set_solver_options( algo="auto", tol=1.0e-8, max_iter=1000, num_steps=NUM_STEPS, dt=TIME_STEP, time_scheme="newmark", newmark_beta=0.25, newmark_gamma=0.5, max_inner_iter=120, inner_tol=1.0e-7, live_viz_interval=2, n_cores=1, device="cpu", ) return session def main() -> None: build_session().compute(mesh_ids=["slider_crank_mesh"]) if __name__ == "__main__": main()