Source code for vivyd.models.rigo_andrianne_denoel

from ..typing import arrf64
from .viv_model import VIVModel
from .generalized_model import GeneralizedModel

from dataclasses import dataclass
from numpy import array, pi




@dataclass(kw_only=True, slots=True)
class RigoAndrianneDenoelModel(VIVModel):
    """
    Rigo, Andrianne & Denoël (2022) wake-oscillator model for the evolution of
    the wake alone, without coupling to the structure :cite:`rigo2022`. The
    equation of the wake is written as
    
    .. math::
        \\ddot{q} - \\left(\\alpha q^2 + \\frac{\\gamma}{\\omega_s^2} \\dot{q}^2 + \\delta \\right) \\omega_s \\dot{q} + \\omega_s^2 q = 0
        
    with
    
    .. math::
        \\omega_s = 2 \\pi \\text{St} \\frac{U_\\infty}{D}.
    
    Parameters
    ----------
    diam : float, optional
        Diameter of the cylinder, :math:`D`. Default is 0.18 m.
    St : float, optional
        Strouhal number, :math:`\\text{St}`. Default is 0.18.
    u_inf : float, optional
        Free-stream velocity, :math:`U_\\infty`. Default is 1.0 m/s.
    alpha : float, optional
        Van der Pol wake nonlinearity parameter, :math:`\\alpha`. Default is -0.085.
    gamma : float, optional
        Rayleigh wake nonlinearity parameter, :math:`\\gamma`. Default is 0.009.
    delta : float, optional
        Constant wake damping parameter, :math:`\\delta`. Default is 0.057.
    """
    state_size = 2
    
    diam : float = 0.18
    """Diameter of the cylinder, :math:`D`."""

    St   : float = 0.18
    """Strouhal number, :math:`\\text{St}`."""

    u_inf: float = 1.0
    """Free-stream velocity, :math:`U_\\infty`."""

    alpha: float = - 0.085
    """Van der Pol wake nonlinearity parameter, :math:`\\alpha`."""

    gamma: float = 0.009
    """Rayleigh wake nonlinearity parameter, :math:`\\gamma`."""

    delta: float = 0.057
    """Constant wake damping parameter, :math:`\\delta`."""
    
    @property
    def St_star(self) -> float:
        """Angular Strouhal number, :math:`\\text{St}^* = 2 \\pi \\text{St}`."""
        
        return 2.0 * pi * self.St
    
    @property
    def ws(self) -> float:
        """Strouhal's vortex shedding pulsation, :math:`\\omega_s = 2 \\pi \\text{St} \\frac{U_\\infty}{D}`."""
        
        return self.St_star * self.u_inf / self.diam
    


    def rhs(self, t: float, state: arrf64, **kwargs) -> arrf64:
        """
        Compute the right-hand side of the system of ODEs.
        
        Args:
            t: Time.
            state: State vector indexed [q_dot, q].

        Returns:
            arrf64: Derivatives of the state vector indexed [dq_dot, dq].
        """
        
        self.validate_state(state)

        q_dot, q = state
        
        dq_dot = (self.alpha * q**2 + self.gamma / self.ws**2 * q_dot**2 + self.delta) * self.ws * q_dot - self.ws**2 * q
        
        dq = q_dot
        
        return array([dq_dot, dq])

    
    @property
    def generalized_params(self) -> dict[str, float]:
        """
        A dictionary containing the parameters of the model in a format
        compatible with the constructor of :class:`GeneralizedModel`.
        """
        return {
            "gamma": self.delta*self.ws,
            "a"    : -self.alpha/self.delta,
            "b"    : -self.gamma/(self.delta*self.ws**2),
            "kappa": self.ws**2
        }
        


    def to_generalized(self) -> GeneralizedModel:
        """
        Returns
        -------
        GeneralizedModel
            A generalized model instance equivalent to the present model.
        """
        return GeneralizedModel(**self.generalized_params)

        


    def to_taichi(self) -> GeneralizedModel:
        """
        Returns
        -------
        GeneralizedModel
            A Taichi-compatible generalized model instance equivalent to the
            present model.
        """
        
        return self.to_generalized()

        


    def get_lift_coefficient(self, state: arrf64, Cl0: float) -> float:
        """
        Compute the lift coefficient from the state vector.
        
        Args:
            state: State vector indexed [q_dot, q].
            Cl0: Lift coefficient at zero angle of attack.
            
        Returns:
            float: Lift coefficient.
        """

        q = state[1]
        return q * Cl0 / 2.0