Dynamic validation example (#111)

examples/DynamicCantileverCase/analytical_dynamic_cantilever.py

@@ -0,0 +1,136 @@
+import numpy as np
+
+
+class AnalyticalDynamicCantilever:
+    """
+    This class computes the analytical solution to a cantilever's dynamic response
+    to an initial velocity profile based on the Euler-Bernoulli beam theory.
+
+    Given the generic dynamic beam equation and boundary conditions imposed on a cantilever,
+    nontrivial solutions exist only when
+
+        cosh(beta * base_length) * cos(beta * base_length) + 1 = 0
+
+    where beta can be used to determine the natural frequencies of the cantilever beam:
+
+        omega = beta^2 * sqrt(E * I / mu)
+
+    The first four roots to beta are given as
+
+        betas = [0.596864*pi, 1.49418*pi, 2.50025*pi, 3.49999*pi] / base_length
+
+    The class solves for the analytical solution at a single natural frequency by computing
+    the non-parametrized mode shapes as well as a mode parameter. The parameter is determined
+    by the initial condition of the rod, namely, the end_velocity input. The free vibration
+    equation can then be applied to determine beam deflection at a given time at any position
+    of the rod:
+
+        w(x, t) = Re[mode_shape * exp(-1j * omega * t)]
+
+    For details, refer to
+    https://en.wikipedia.org/wiki/Euler%E2%80%93Bernoulli_beam_theory#Dynamic_beam_equation
+    or
+    Han et al (1998), Dynamics of Transversely Vibrating Beams
+    Using Four Engineering Theories
+
+        Attributes
+        ----------
+        base_length: float
+            Total length of the rod
+        base_area: float
+            Cross-sectional area of the rod
+        moment_of_inertia: float
+            Second moment of area of the rod's cross-section
+        young's_modulus: float
+            Young's modulus of the rod
+        density: float
+            Density of the rod
+        mode: int
+            Index of the first 'mode' th natural frequency.
+            Up to the first four modes are supported.
+        end_velocity: float
+            Initial oscillatory velocity at the end of the rod
+    """
+
+    def __init__(
+        self,
+        base_length,
+        base_area,
+        moment_of_inertia,
+        youngs_modulus,
+        density,
+        mode,
+        end_velocity=0.0,
+    ):
+        self.base_length = base_length
+        self.end_velocity = end_velocity
+
+        if not (isinstance(mode, int) and mode >= 0 and mode < 5):
+            raise ValueError(
+                "Unsupported mode value, please provide a integer value from 0-4"
+            )
+
+        self.mode = mode
+
+        betas = np.array([0.596864, 1.49418, 2.50025, 3.49999]) * np.pi / base_length
+        self.beta = betas[mode]
+        self.omega = (self.beta ** 2) * np.sqrt(
+            youngs_modulus
+            * moment_of_inertia
+            / (density * base_area * base_length ** 4)
+        )
+
+        nonparametrized_mode_at_end = self._compute_nonparametrized_mode(
+            base_length, base_length, self.beta
+        )
+        self.mode_param = end_velocity / (
+            -1j * self.omega * nonparametrized_mode_at_end
+        )
+
+    def get_initial_velocity_profile(self, positions):
+        initial_velocities = []
+
+        for x in positions:
+            initial_velocities.append(
+                np.real(
+                    (-1j * self.omega * self.mode_param)
+                    * self._compute_nonparametrized_mode(
+                        self.base_length * x, self.base_length, self.beta
+                    )
+                )
+            )
+        return np.array(initial_velocities)
+
+    def get_time_dependent_positions(self, position, time_array):
+        time_dependent_positions = (
+            self.mode_param
+            * self._compute_nonparametrized_mode(
+                self.base_length * position, self.base_length, self.beta
+            )
+            * np.exp(-1j * self.omega * time_array)
+        )
+        return np.real(time_dependent_positions)
+
+    def get_time_dependent_velocities(self, position, time_array):
+        time_dependent_velocities = (
+            (-1j * self.omega * self.mode_param)
+            * self._compute_nonparametrized_mode(
+                self.base_length * position, self.base_length, self.beta
+            )
+            * np.exp(-1j * self.omega * time_array)
+        )
+        return np.real(time_dependent_velocities)
+
+    def get_omega(self):
+        return self.omega
+
+    def get_amplitude(self):
+        return self.end_velocity / self.omega
+
+    @staticmethod
+    def _compute_nonparametrized_mode(x, base_length, beta):
+        a = np.cosh(beta * x) - np.cos(beta * x)
+        b = np.cos(beta * base_length) + np.cosh(beta * base_length)
+        c = np.sin(beta * base_length) + np.sinh(beta * base_length)
+        d = np.sin(beta * x) - np.sinh(beta * x)
+        return a + b / c * d

examples/DynamicCantileverCase/dynamic_cantilever.py

@@ -0,0 +1,170 @@
+import numpy as np
+from scipy.fft import fft, fftfreq
+from scipy.signal import find_peaks
+from elastica import *
+from analytical_dynamic_cantilever import AnalyticalDynamicCantilever
+
+
+def simulate_dynamic_cantilever_with(
+    density=2000.0,
+    n_elem=100,
+    final_time=300.0,
+    mode=0,
+    rendering_fps=30.0,  # For visualization
+):
+    """
+    This function completes a dynamic cantilever simulation with the given parameters.
+
+    Parameters
+    ----------
+    density: float
+        Density of the rod
+    n_elem: int
+        The number of elements of the rod
+    final_time: float
+        Total simulation time. The timestep is determined by final_time / n_steps.
+    mode: int
+        Index of the first 'mode' th natural frequency.
+        Up to the first four modes are supported.
+    rendering_fps: float
+        Frames per second for video plotting.
+        The call back step-skip is also determined by rendering_fps.
+
+    Returns
+    -------
+    dict of {str : int}
+        A collection of parameters for post-processing.
+
+    """
+
+    class DynamicCantileverSimulator(BaseSystemCollection, Constraints, CallBacks):
+        pass
+
+    cantilever_sim = DynamicCantileverSimulator()
+
+    # Add test parameters
+    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
+    base_radius = 0.02
+    base_area = np.pi * base_radius ** 2
+    youngs_modulus = 1e5
+
+    moment_of_inertia = np.pi / 4 * base_radius ** 4
+
+    dl = base_length / n_elem
+    dt = dl * 0.05
+    step_skips = int(1.0 / (rendering_fps * dt))
+
+    # Add Cosserat rod
+    cantilever_rod = CosseratRod.straight_rod(
+        n_elem,
+        start,
+        direction,
+        normal,
+        base_length,
+        base_radius,
+        density,
+        0.0,
+        youngs_modulus,
+    )
+
+    # Add constraints
+    cantilever_sim.append(cantilever_rod)
+    cantilever_sim.constrain(cantilever_rod).using(
+        OneEndFixedBC, constrained_position_idx=(0,), constrained_director_idx=(0,)
+    )
+
+    end_velocity = 0.005
+    analytical_cantilever_soln = AnalyticalDynamicCantilever(
+        base_length,
+        base_area,
+        moment_of_inertia,
+        youngs_modulus,
+        density,
+        mode=mode,
+        end_velocity=end_velocity,
+    )
+
+    initial_velocity = analytical_cantilever_soln.get_initial_velocity_profile(
+        cantilever_rod.position_collection[0, :]
+    )
+    cantilever_rod.velocity_collection[2, :] = initial_velocity
+
+    # Add call backs
+    class CantileverCallBack(CallBackBaseClass):
+        def __init__(self, step_skip: int, callback_params: dict):
+            CallBackBaseClass.__init__(self)
+            self.every = step_skip
+            self.callback_params = callback_params
+
+        def make_callback(self, system, time, current_step: int):
+
+            if current_step % self.every == 0:
+
+                self.callback_params["time"].append(time)
+                self.callback_params["position"].append(
+                    system.position_collection.copy()
+                )
+                self.callback_params["deflection"].append(
+                    system.position_collection[2, -1].copy()
+                )
+                return
+
+    recorded_history = defaultdict(list)
+    cantilever_sim.collect_diagnostics(cantilever_rod).using(
+        CantileverCallBack, step_skip=step_skips, callback_params=recorded_history
+    )
+    cantilever_sim.finalize()
+
+    total_steps = int(final_time / dt)
+    print(f"Total steps: {total_steps}")
+
+    timestepper = PositionVerlet()
+
+    integrate(
+        timestepper,
+        cantilever_sim,
+        final_time,
+        total_steps,
+    )
+
+    # FFT
+    amplitudes = np.abs(fft(recorded_history["deflection"]))
+    fft_length = len(amplitudes)
+    amplitudes = amplitudes * 2 / fft_length
+    omegas = fftfreq(fft_length, dt * step_skips) * 2 * np.pi  # [rad/s]
+
+    try:
+        peaks, _ = find_peaks(amplitudes)
+        peak = peaks[np.argmax(amplitudes[peaks])]
+
+        simulated_frequency = omegas[peak]
+        theoretical_frequency = analytical_cantilever_soln.get_omega()
+
+        simulated_amplitude = max(recorded_history["deflection"])
+        theoretical_amplitude = analytical_cantilever_soln.get_amplitude()
+
+        print(
+            f"Theoretical frequency: {theoretical_frequency} rad/s \n"
+            f"Simulated frequency: {simulated_frequency} rad/s \n"
+            f"Theoretical amplitude: {theoretical_amplitude} m \n"
+            f"Simulated amplitude: {simulated_amplitude} m"
+        )
+
+        return {
+            "rod": cantilever_rod,
+            "recorded_history": recorded_history,
+            "fft_frequencies": omegas,
+            "fft_amplitudes": amplitudes,
+            "analytical_cantilever_soln": analytical_cantilever_soln,
+            "peak": peak,
+            "simulated_frequency": simulated_frequency,
+            "theoretical_frequency": theoretical_frequency,
+            "simulated_amplitude": simulated_amplitude,
+            "theoretical_amplitude": theoretical_amplitude,
+        }
+
+    except RuntimeError:
+        print("No peaks detected: change input parameters.")

examples/DynamicCantileverCase/dynamic_cantilever_phase_space.py

@@ -0,0 +1,40 @@
+__doc__ = """ Validating phase space of dynamic cantilever beam analytical_cantilever_soln with  respect to varying densities.
+The theoretical dynamic response is obtained via Euler-Bernoulli beam theory."""
+
+from dynamic_cantilever_post_processing import plot_phase_space_with
+from dynamic_cantilever import simulate_dynamic_cantilever_with
+
+if __name__ == "__main__":
+    import multiprocessing as mp
+
+    PLOT_FIGURE = True
+    SAVE_FIGURE = False
+
+    # density = 500, 1000, 2000, 3000, 4000, 5000 kg/m3
+    densities = [500]
+    densities.extend(range(1000, 6000, 1000))
+
+    with mp.Pool(mp.cpu_count()) as pool:
+        results = pool.map(simulate_dynamic_cantilever_with, densities)
+
+    theory_frequency = []
+    sim_frequency = []
+    theory_amplitude = []
+    sim_amplitude = []
+
+    for result in results:
+        theory_frequency.append(result["theoretical_frequency"])
+        sim_frequency.append(result["simulated_frequency"])
+        theory_amplitude.append(result["theoretical_amplitude"])
+        sim_amplitude.append(result["simulated_amplitude"])
+
+    # Plot frequencies and amplitudes vs densities
+    plot_phase_space_with(
+        densities,
+        theory_frequency,
+        sim_frequency,
+        theory_amplitude,
+        sim_amplitude,
+        PLOT_FIGURE=PLOT_FIGURE,
+        SAVE_FIGURE=SAVE_FIGURE,
+    )