.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "_gallery/AxialStretchingCase/run_axial_stretching.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr__gallery_AxialStretchingCase_run_axial_stretching.py: Axial Stretching ================ This case tests the axial stretching of a rod. The expected behavior is supposed to be like a spring-gravity motion, but with a rod. A rod is fixed at one end and a force is applied at the other end. The rod stretches and the displacement of the tip is compared with the analytical solution. .. GENERATED FROM PYTHON SOURCE LINES 11-20 .. code-block:: Python # isort:skip_file import numpy as np from collections import defaultdict from matplotlib import pyplot as plt import elastica as ea .. GENERATED FROM PYTHON SOURCE LINES 21-25 Simulation Setup ---------------- We define a simulator class that inherits from the necessary mixins. This makes constraints, forces, and damping available to the system. .. GENERATED FROM PYTHON SOURCE LINES 25-36 .. code-block:: Python class StretchingBeamSimulator( ea.BaseSystemCollection, ea.Constraints, ea.Forcing, ea.Damping, ea.CallBacks ): pass stretch_sim = StretchingBeamSimulator() final_time = 200.0 .. GENERATED FROM PYTHON SOURCE LINES 37-44 Rod Setup --------- Next, we set up the test parameters for the simulating rods. This includes the number of elements, the start position, direction, normal, length, radius, density, and Young's modulus of the rod. For this case, we have fixed boundary condition at one end, and we apply external force at the other end. .. GENERATED FROM PYTHON SOURCE LINES 44-83 .. code-block:: Python # setting up test params n_elem = 19 start = np.zeros((3,)) direction = np.array([1.0, 0.0, 0.0]) normal = np.array([0.0, 1.0, 0.0]) base_length = 1.0 base_radius = 0.025 base_area = np.pi * base_radius**2 density = 1000 youngs_modulus = 1e4 # For shear modulus of 1e4, nu is 99! poisson_ratio = 0.5 shear_modulus = youngs_modulus / (poisson_ratio + 1.0) stretchable_rod = ea.CosseratRod.straight_rod( n_elem, start, direction, normal, base_length, base_radius, density, youngs_modulus=youngs_modulus, shear_modulus=shear_modulus, ) stretch_sim.append(stretchable_rod) stretch_sim.constrain(stretchable_rod).using( ea.OneEndFixedBC, constrained_position_idx=(0,), constrained_director_idx=(0,) ) end_force_x = 1.0 end_force = np.array([end_force_x, 0.0, 0.0]) stretch_sim.add_forcing_to(stretchable_rod).using( ea.EndpointForces, 0.0 * end_force, end_force, ramp_up_time=1e-2 ) .. GENERATED FROM PYTHON SOURCE LINES 84-86 Damping is added to the system to help it reach a steady state. We use an `AnalyticalLinearDamper` to add damping to the rod. .. GENERATED FROM PYTHON SOURCE LINES 86-98 .. code-block:: Python # add damping dl = base_length / n_elem dt = 0.1 * dl damping_constant = 0.1 stretch_sim.dampen(stretchable_rod).using( ea.AnalyticalLinearDamper, damping_constant=damping_constant, time_step=dt, ) .. GENERATED FROM PYTHON SOURCE LINES 99-103 Callbacks --------- A callback object is passed to the simulator to record states of the rod during the simulation. This is useful for post-processing the results. .. GENERATED FROM PYTHON SOURCE LINES 103-138 .. code-block:: Python # Add call backs class AxialStretchingCallBack(ea.CallBackBaseClass): """ Tracks the velocity norms of the rod """ def __init__(self, step_skip: int, callback_params: dict) -> None: super().__init__() self.every = step_skip self.callback_params = callback_params def make_callback( self, system: ea.typing.RodType, time: np.float64, current_step: int ) -> None: if current_step % self.every == 0: self.callback_params["time"].append(time) # Collect only x self.callback_params["position"].append( system.position_collection[0, -1].copy() ) self.callback_params["velocity_norms"].append( np.linalg.norm(system.velocity_collection.copy()) ) return recorded_history: dict[str, list] = defaultdict(list) stretch_sim.collect_diagnostics(stretchable_rod).using( AxialStretchingCallBack, step_skip=200, callback_params=recorded_history ) .. GENERATED FROM PYTHON SOURCE LINES 139-143 Finalize and Run ---------------- We finalize the simulator and create the time-stepper. The `PositionVerlet` time-stepper is used to integrate the system. .. GENERATED FROM PYTHON SOURCE LINES 143-155 .. code-block:: Python stretch_sim.finalize() timestepper: ea.typing.StepperProtocol = ea.PositionVerlet() # timestepper = PEFRL() total_steps = int(final_time / dt) print("Total steps", total_steps) dt = final_time / total_steps time = 0.0 for i in range(total_steps): time = timestepper.step(stretch_sim, time, dt) .. rst-class:: sphx-glr-script-out .. code-block:: none Total steps 38000 .. GENERATED FROM PYTHON SOURCE LINES 156-161 Post-Processing --------------- Finally, we plot the results and compare them with the analytical solution. The analytical solution is calculated using the first-order theory with both the base length and the modified length. .. GENERATED FROM PYTHON SOURCE LINES 161-177 .. code-block:: Python # First-order theory with base-length expected_tip_disp = end_force_x * base_length / base_area / youngs_modulus # First-order theory with modified-length, gives better estimates expected_tip_disp_improved = ( end_force_x * base_length / (base_area * youngs_modulus - end_force_x) ) fig = plt.figure(figsize=(10, 8), frameon=True, dpi=150) ax = fig.add_subplot(111) ax.plot(recorded_history["time"], recorded_history["position"], lw=2.0) ax.hlines(base_length + expected_tip_disp, 0.0, final_time, "k", "dashdot", lw=1.0) ax.hlines( base_length + expected_tip_disp_improved, 0.0, final_time, "k", "dashed", lw=2.0 ) plt.show() .. image-sg:: /_gallery/AxialStretchingCase/images/sphx_glr_run_axial_stretching_001.png :alt: run axial stretching :srcset: /_gallery/AxialStretchingCase/images/sphx_glr_run_axial_stretching_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 3.160 seconds) .. _sphx_glr_download__gallery_AxialStretchingCase_run_axial_stretching.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: run_axial_stretching.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: run_axial_stretching.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: run_axial_stretching.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_