conspire/constitutive/fluid/viscoplastic/

mod.rs

//! Viscoplastic fluid constitutive models.

#[cfg(feature = "autodiff")]
pub mod autodiff;

#[cfg(test)]
mod test;

use crate::{
    constitutive::{ConstitutiveError, fluid::plastic::Plastic},
    math::{
        Derivative, Differentiable, IDENTITY_22, Intermediate, Quantity, Rank2, Scalar, Tensor,
        TensorArray, TensorRank4, TensorTuple, TensorTupleVec,
    },
    mechanics::{
        DeformationGradientPlastic, DeformationGradientRatePlastic, MandelStressElastic,
        StretchingRatePlastic, StretchingRatePlasticTangent, StretchingRatePlasticTangentYield,
    },
    units::{Dissipation, Rate, Stress},
};

type IntermediateRank4 = TensorRank4<3, Intermediate, Intermediate, Intermediate, Intermediate>;

/// Viscoplastic state variables.
pub type ViscoplasticStateVariables<Y> = TensorTuple<DeformationGradientPlastic, Y>;

/// Viscoplastic state variables history.
pub type ViscoplasticStateVariablesHistory<Y> = TensorTupleVec<DeformationGradientPlastic, Y>;

/// The evolution of the viscoplastic state variables.
pub type ViscoplasticEvolution<Y> = Derivative<ViscoplasticStateVariables<Y>>;

/// The viscoplastic state's Lie-algebra rate `(D_p, Ẏ)`, as consumed by the
/// field time-integration drivers — the plastic stretching rate itself, not the
/// group velocity `Ḟ_p = D_p F_p`, paired with the hardening-variable rate.
pub type ViscoplasticAlgebraRate<Y = Quantity> = TensorTuple<StretchingRatePlastic, Derivative<Y>>;

/// The history of the evolution of the viscoplastic state variables.
pub type ViscoplasticEvolutionHistory<Y> =
    TensorTupleVec<DeformationGradientRatePlastic, Derivative<Y>>;

/// Required methods for viscoplastic fluid constitutive models.
pub trait Viscoplastic<Y>
where
    Self: Plastic,
    Y: Differentiable + Tensor,
{
    /// Returns the initial state of the variables.
    fn initial_state(&self) -> ViscoplasticStateVariables<Y>;
    /// Calculates and returns the plastic evolution.
    ///
    /// ```math
    /// \dot{\mathbf{F}}_\mathrm{p} = \mathbf{D}_\mathrm{p}\cdot\mathbf{F}_\mathrm{p}\quad\text{and}\quad\dot{\varepsilon}_\mathrm{p} = |\mathbf{D}_\mathrm{p}|
    /// ```
    fn plastic_evolution(
        &self,
        mandel_stress: MandelStressElastic,
        state_variables: &ViscoplasticStateVariables<Y>,
    ) -> Result<ViscoplasticEvolution<Y>, ConstitutiveError>;
    /// Calculates and returns the rate of plastic stretching.
    ///
    /// ```math
    /// \mathbf{D}_\mathrm{p}(\mathbf{M}_\mathrm{e}') = d_0\left(\frac{|\mathbf{M}_\mathrm{e}'|}{Y(S)}\right)^{\footnotesize\tfrac{1}{m}}\frac{\mathbf{M}_\mathrm{e}'}{|\mathbf{M}_\mathrm{e}'|}
    /// ```
    fn plastic_stretching_rate(
        &self,
        deviatoric_mandel_stress: MandelStressElastic,
        yield_stress: Quantity<Stress>,
    ) -> Result<StretchingRatePlastic, ConstitutiveError> {
        let magnitude = deviatoric_mandel_stress.norm();
        if magnitude.is_zero() {
            Ok(StretchingRatePlastic::zero())
        } else {
            let reference_flow_rate = self.reference_flow_rate();
            Ok(deviatoric_mandel_stress
                * (reference_flow_rate / magnitude
                    * (magnitude / yield_stress).powf(1.0 / self.rate_sensitivity())))
        }
    }
    /// Calculates and returns the tangent of the plastic stretching rate with
    /// respect to the deviatoric Mandel stress.
    ///
    /// ```math
    /// \frac{\partial D^\mathrm{p}_{ij}}{\partial M'_{kl}} = g\,\delta_{ik}\delta_{jl} + g\left(\frac{1}{m} - 1\right)\frac{M'_{ij}M'_{kl}}{|\mathbf{M}'|^2},
    /// \qquad g = \frac{d_0}{|\mathbf{M}'|}\left(\frac{|\mathbf{M}'|}{Y}\right)^{\footnotesize\tfrac{1}{m}}
    /// ```
    fn plastic_stretching_rate_tangent(
        &self,
        deviatoric_mandel_stress: &MandelStressElastic,
        yield_stress: Quantity<Stress>,
    ) -> Result<StretchingRatePlasticTangent, ConstitutiveError> {
        let magnitude = deviatoric_mandel_stress.norm();
        if magnitude.is_zero() {
            Ok(StretchingRatePlasticTangent::zero())
        } else {
            let rate_sensitivity = self.rate_sensitivity();
            let flow = self.reference_flow_rate() / magnitude
                * (magnitude / yield_stress).powf(1.0 / rate_sensitivity);
            let normal = deviatoric_mandel_stress / magnitude;
            Ok(
                IntermediateRank4::dyad_ik_jl(&IDENTITY_22, &IDENTITY_22) * flow
                    + IntermediateRank4::dyad_ij_kl(&normal, &normal)
                        * (flow * (1.0 / rate_sensitivity - 1.0)),
            )
        }
    }
    /// Calculates and returns the tangent of the plastic stretching rate with
    /// respect to the yield stress.
    ///
    /// ```math
    /// \frac{\partial\mathbf{D}^\mathrm{p}}{\partial Y} = -\frac{\mathbf{D}^\mathrm{p}}{mY}
    /// ```
    fn plastic_stretching_rate_tangent_yield(
        &self,
        deviatoric_mandel_stress: MandelStressElastic,
        yield_stress: Quantity<Stress>,
    ) -> Result<StretchingRatePlasticTangentYield, ConstitutiveError> {
        Ok(
            self.plastic_stretching_rate(deviatoric_mandel_stress, yield_stress)? / yield_stress
                * (-1.0 / self.rate_sensitivity()),
        )
    }
    /// Calculates and returns the dissipation potential.
    ///
    /// ```math
    /// \phi(\mathbf{D}_\mathrm{p}) = \frac{Yd_0}{1+m}\left(\frac{|\mathbf{D}_\mathrm{p}|}{d_0}\right)^{1+m}
    /// ```
    fn dissipation_potential(
        &self,
        plastic_stretching_rate: StretchingRatePlastic,
        yield_stress: Quantity<Stress>,
    ) -> Result<Quantity<Dissipation>, ConstitutiveError> {
        let rate_sensitivity = self.rate_sensitivity();
        let reference_flow_rate = self.reference_flow_rate();
        Ok(
            reference_flow_rate * yield_stress / (1.0 + rate_sensitivity)
                * (plastic_stretching_rate.norm() / reference_flow_rate)
                    .powf(1.0 + rate_sensitivity),
        )
    }
    /// Calculates and returns the dual dissipation potential.
    ///
    /// ```math
    /// \phi^*(\mathbf{M}_\mathrm{e}') = \frac{Yd_0m}{1+m}\left(\frac{|\mathbf{M}_\mathrm{e}'|}{Y(S)}\right)^{\footnotesize\tfrac{1+m}{m}}
    /// ```
    fn dual_dissipation_potential(
        &self,
        deviatoric_mandel_stress: MandelStressElastic,
        yield_stress: Quantity<Stress>,
    ) -> Result<Quantity<Dissipation>, ConstitutiveError> {
        let rate_sensitivity = self.rate_sensitivity();
        Ok(self.reference_flow_rate()
            * yield_stress
            * (rate_sensitivity / (1.0 + rate_sensitivity))
            * (deviatoric_mandel_stress.norm() / yield_stress)
                .powf((1.0 + rate_sensitivity) / rate_sensitivity))
    }
    /// Returns the rate_sensitivity parameter.
    fn rate_sensitivity(&self) -> Scalar;
    /// Returns the reference flow rate.
    fn reference_flow_rate(&self) -> Quantity<Rate>;
}

/// The viscoplastic flow model.
#[derive(Clone, Debug)]
pub struct ViscoplasticFlow {
    /// The initial yield stress $`Y_0`$.
    pub yield_stress: Quantity<Stress>,
    /// The isotropic hardening slope $`H`$.
    pub hardening_slope: Quantity<Stress>,
    /// The rate sensitivity parameter $`m`$.
    pub rate_sensitivity: Scalar,
    /// The reference flow rate $`d_0`$.
    pub reference_flow_rate: Quantity<Rate>,
}

impl Plastic for ViscoplasticFlow {
    fn initial_yield_stress(&self) -> Quantity<Stress> {
        self.yield_stress
    }
    fn hardening_slope(&self) -> Quantity<Stress> {
        self.hardening_slope
    }
}

impl Viscoplastic<Quantity> for ViscoplasticFlow {
    fn initial_state(&self) -> ViscoplasticStateVariables<Quantity> {
        (DeformationGradientPlastic::identity(), Quantity::default()).into()
    }
    fn plastic_evolution(
        &self,
        mandel_stress: MandelStressElastic,
        state_variables: &ViscoplasticStateVariables<Quantity>,
    ) -> Result<ViscoplasticEvolution<Quantity>, ConstitutiveError> {
        default_plastic_evolution(self, mandel_stress, state_variables)
    }
    fn rate_sensitivity(&self) -> Scalar {
        self.rate_sensitivity
    }
    fn reference_flow_rate(&self) -> Quantity<Rate> {
        self.reference_flow_rate
    }
}

pub fn default_plastic_evolution<C>(
    model: &C,
    mandel_stress: MandelStressElastic,
    state_variables: &ViscoplasticStateVariables<Quantity>,
) -> Result<ViscoplasticEvolution<Quantity>, ConstitutiveError>
where
    C: Viscoplastic<Quantity>,
{
    let (deformation_gradient_p, &equivalent_plastic_strain) = state_variables.into();
    let plastic_stretching_rate = model.plastic_stretching_rate(
        mandel_stress.deviatoric(),
        model.yield_stress(equivalent_plastic_strain)?,
    )?;
    let equivalent_plastic_strain_rate = plastic_stretching_rate.norm();
    Ok((
        plastic_stretching_rate * deformation_gradient_p,
        equivalent_plastic_strain_rate,
    )
        .into())
}