Implementation of plastic dilation

Author: YiminJin

doc/modules/changes/20250616_yiminjin

@@ -0,0 +1,11 @@
+Changed: The prescribed dilation is now splitted into two parts: a
+pressure-dependent part that is moved to the left-hand side of the mass
+conservation equation, and a pressure-independent part that stays on the
+right-hand side. In addition, the effect of dilation on the momentum
+conservation equation for incompressible models is no longer represented by
+the dilation term on the right-hand side, but by the volumetric strain rate
+term on the left-hand side. These changes enable the Drucker-Prager model
+to produce associated plastic flows, in which the shear bands develop along
+the Coulomb angle.
+<br>
+(Yimin Jin, 2025/06/16)

include/aspect/material_model/interface.h

@@ -1204,15 +1204,27 @@ namespace aspect
      * Stokes solution.
      *
      * This is typically used in a MaterialModel to add dilation when plastic
-     * failure occurs as motivated by ChoiPeterson2015. If this output
-     * (denoted by R below) is present and enable_prescribed_dilation==true
-     * the following terms will be assembled:
-     *
-     * 1) $\int - (R,q)$ to the conservation of mass equation, creating
-     *    $-(div u,q) = -(R,q)$.
-     * 2) $\int - 2.0 / 3.0 * eta * (R, div v)$ to the RHS of the momentum
-     *    equation (if the model is incompressible), otherwise this term is
-     *    already present on the left side.
+     * failure occurs, but can be used for other forms of dilation as well.
+     * When plastic dilation is included, a term
+     * $\bar\alpha\gamma$ should be added to the right-hand side of the
+     * mass conservation equation, where $\bar\alpha$ is the negative
+     * derivative of plastic potential with respect to the pressure
+     * ($\sin\psi$ in 2D case), and $\gamma$ is the plastic multiplier.
+     * The plastic multiplier is given by
+     * $\gamma = (\tau_{II} - \alpha p - k) / \eta^{ve}$,
+     * where $\tau_{II}$ is the second invariant of the deviatoric stress,
+     * $\alpha$ is the negative derivative of yield function with respect to
+     * the pressure ($\sin\phi$ in 2D case), $k$ is cohesion, and $\eta^{ve}$,
+     * is the pre-yielding viscosity. When the Picard method or Defect Correction
+     * Method is applied, the term $\bar\alpha\gamma$ should be split into two
+     * terms:
+     * $\bar\alpha\gamma = \bar\alpha\alpha p / \eta^{ve} +
+     * \bar\alpha(\tau_{II} - k) / \eta^{ve}$,
+     * the former of which should be moved to the left-hand side in order to
+     * guarantee the stability of the nonlinear solver. Therefore, this output
+     * provides two terms: dilation_lhs_term corresponds to
+     * $\bar\alpha\alpha / \eta^{ve}$ (p is replaced by the shape function),
+     * and dilation_rhs_term corresponds to $(\tau_{II} - k) / \eta^{ve}$.
      */
     template <int dim>
     class PrescribedPlasticDilation : public NamedAdditionalMaterialOutputs<dim>
@@ -1229,10 +1241,16 @@ namespace aspect
         std::vector<double> get_nth_output(const unsigned int idx) const override;
 
         /**
-         * A scalar value per evaluation point that specifies the prescribed dilation
-         * in that point.
+         * A scalar value per evaluation point corresponding to the LHS term
+         * due to plastic dilation.
+         */
+        std::vector<double> dilation_lhs_term;
+
+        /**
+         * A scalar value per evaluation point corresponding to the RHS term
+         * due to plastic dilation.
          */
-        std::vector<double> dilation;
+        std::vector<double> dilation_rhs_term;
     };
 
 

include/aspect/material_model/rheology/drucker_prager.h

@@ -43,6 +43,11 @@ namespace aspect
          */
         double angle_internal_friction;
 
+        /**
+         * Dilation angle (psi) for the current composition and phase
+         */
+        double angle_dilation;
+
         /**
          * Cohesion for the current composition and phase
          */
@@ -170,6 +175,20 @@ namespace aspect
           compute_derivative (const double angle_internal_friction,
                               const double effective_strain_rate) const;
 
+          /**
+           * Compute the LHS and RHS dilation terms for the Stokes system.
+           * LHS: $\bar\alpha\alpha / \eta^{ve}$;
+           * RHS: $\bar(2\eta^{ve}\varepsilon^{eff} - k) / \eta^{ve}$.
+           * Here $\alpha$ and $\bar\alpha$ correspond to the friction angle
+           * and the dilation angle, respectively, $k$ is cohesion,
+           * $\eta^{ve}$ is the non-yielding viscosity, and $\varepsilon^{eff}$
+           * is the effective viscosity.
+           */
+          std::pair<double,double>
+          compute_dilation_terms_for_stokes_system(const DruckerPragerParameters &drucker_prager_parameters,
+                                                   const double non_yielding_viscosity,
+                                                   const double effective_strain_rate) const;
+
         private:
 
           /**
@@ -183,6 +202,17 @@ namespace aspect
            */
           std::vector<double> angles_internal_friction;
 
+          /**
+           * The plastic potential of the Drucker-Prager model is given by
+           * tau_II - 6.0 * pressure * sin_psi / (sqrt(3.0) * (3.0 + sin_psi)) in 3D or
+           * tau_II - pressure * sin_psi in 2D, where tau_II is the second invariant
+           * of the deviatoric stress.
+           * Psi is an angle of dilation, that is input by the user in degrees,
+           * but stored as radians. Note that the dilation angle must not be larger
+           * than the internal friction angle.
+           */
+          std::vector<double> angles_dilation;
+
           /**
            * The cohesion is provided and stored in Pa.
            */

include/aspect/material_model/rheology/visco_plastic.h

@@ -102,6 +102,20 @@ namespace aspect
        * All the drucker prager plasticity parameters.
        */
       std::vector<Rheology::DruckerPragerParameters> drucker_prager_parameters;
+
+      /**
+       * The LHS term corresponding to plastic dilation in the
+       * Stokes system. For details, see the comments of
+       * MaterialModel::PrescribedDilation::dilation_lhs_term.
+       */
+      std::vector<double> dilation_lhs_terms;
+
+      /**
+       * The RHS term corresponding to plastic dilation in the
+       * Stokes system. For details, see the comments of
+       * MaterialModel::PrescribedDilation::dilation_rhs_term.
+       */
+      std::vector<double> dilation_rhs_terms;
     };
 
     namespace Rheology
@@ -143,7 +157,7 @@ namespace aspect
            */
           void compute_viscosity_derivatives(const unsigned int point_index,
                                              const std::vector<double> &volume_fractions,
-                                             const std::vector<double> &composition_viscosities,
+                                             const IsostrainViscosities &isostrain_values,
                                              const MaterialModel::MaterialModelInputs<dim> &in,
                                              MaterialModel::MaterialModelOutputs<dim> &out,
                                              const std::vector<double> &phase_function_values = std::vector<double>(),

include/aspect/newton.h

@@ -61,6 +61,19 @@ namespace aspect
          * derivatives when material averaging is applied.
          */
         std::vector<double> viscosity_derivative_averaging_weights;
+
+        /**
+         * The derivatives of the plastic dilation with respect to pressure.
+         * Note that this term does not include $\bar\alpha\alpha / \eta^{ve}$,
+         * which is always on the left-hand side no matter if the Newton method
+         * is applied or not.
+         */
+        std::vector<double> dilation_derivative_wrt_pressure;
+
+        /**
+         * The derivatives of the plastic dilation with respect to strain rate.
+         */
+        std::vector<SymmetricTensor<2,dim>> dilation_derivative_wrt_strain_rate;
     };
   }
 

include/aspect/simulator/assemblers/stokes.h

@@ -40,6 +40,8 @@ namespace aspect
         void
         execute(internal::Assembly::Scratch::ScratchBase<dim>   &scratch_base,
                 internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const override;
+
+        void create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &outputs) const override;
     };
 
     /**

include/aspect/simulator/solver/matrix_free_operators.h

@@ -107,6 +107,11 @@ namespace aspect
        */
       bool enable_newton_derivatives;
 
+      /**
+       * If true, plastic dilation terms are part of the operator.
+       */
+      bool enable_prescribed_dilation;
+
       /**
        * Symmetrize the Newton system when it's true (i.e., the
        * stabilization is symmetric or SPD).
@@ -149,6 +154,30 @@ namespace aspect
       Table<2, SymmetricTensor<2, dim, VectorizedArray<number>>>
       newton_factor_wrt_strain_rate_table;
 
+      /**
+       * Table which stores the dilation LHS terms. For the definition of this
+       * term, see the comments of MaterialModel::PrescribedPlasticDilation::
+       * dilation_lhs_term.
+       */
+      Table<2, VectorizedArray<number>> dilation_lhs_term_table;
+
+      /**
+       * Table which stores the product of the plastic dilation derivative
+       * with respect to pressure and the Newton derivative scaling factor.
+       * Note that the plastic dilation derivative does not include the term
+       * $\bar\alpha\alpha / \eta^{ve}$, which is always on the left-hand
+       * side no matter if the Newton method is applied.
+       */
+      Table<2, VectorizedArray<number>>
+      dilation_derivative_wrt_pressure_table;
+
+      /**
+       * Table which stores the product of the plastic dilation derivative
+       * with respect to strain-rate and the Newton derivative scaling factor.
+       */
+      Table<2, SymmetricTensor<2, dim, VectorizedArray<number>>>
+      dilation_derivative_wrt_strain_rate_table;
+
       /**
        * Table which stores the product of the pressure perturbation
        * and the normalized gravity. The size is n_face_boundary * n_face_q_points,

source/material_model/interface.cc

@@ -998,20 +998,45 @@ namespace aspect
 
 
 
+    namespace
+    {
+      std::vector<std::string> make_prescribed_dilation_outputs_names()
+      {
+        std::vector<std::string> names;
+        names.emplace_back("dilation_lhs_term");
+        names.emplace_back("dilation_rhs_term");
+        return names;
+      }
+    }
+
+
+
     template <int dim>
     PrescribedPlasticDilation<dim>::PrescribedPlasticDilation (const unsigned int n_points)
-      : NamedAdditionalMaterialOutputs<dim>(std::vector<std::string>(1, "prescribed_dilation")),
-        dilation(n_points, numbers::signaling_nan<double>())
+      : NamedAdditionalMaterialOutputs<dim>(make_prescribed_dilation_outputs_names()),
+        dilation_lhs_term(n_points, numbers::signaling_nan<double>()),
+        dilation_rhs_term(n_points, numbers::signaling_nan<double>())
     {}
 
 
 
     template <int dim>
     std::vector<double> PrescribedPlasticDilation<dim>::get_nth_output(const unsigned int idx) const
     {
-      (void)idx;
-      Assert(idx==0, ExcInternalError());
-      return dilation;
+      AssertIndexRange (idx, 2);
+      switch (idx)
+        {
+          case 0:
+            return dilation_lhs_term;
+
+          case 1:
+            return dilation_rhs_term;
+
+          default:
+            AssertThrow(false, ExcInternalError());
+        }
+      // we will never get here, so just return something
+      return std::vector<double>();
     }
 
 

source/material_model/rheology/drucker_prager.cc

@@ -34,7 +34,8 @@ namespace aspect
     {
       DruckerPragerParameters::DruckerPragerParameters()
         : angle_internal_friction (numbers::signaling_nan<double>()),
-          cohesion  (numbers::signaling_nan<double>()),
+          angle_dilation (numbers::signaling_nan<double>()),
+          cohesion (numbers::signaling_nan<double>()),
           max_yield_stress (numbers::signaling_nan<double>())
       {}
 
@@ -58,6 +59,7 @@ namespace aspect
           {
             // no phases
             drucker_prager_parameters.angle_internal_friction = angles_internal_friction[composition];
+            drucker_prager_parameters.angle_dilation = angles_dilation[composition];
             drucker_prager_parameters.cohesion = cohesions[composition];
             drucker_prager_parameters.max_yield_stress = max_yield_stresses[composition];
           }
@@ -66,6 +68,8 @@ namespace aspect
             // Average among phases
             drucker_prager_parameters.angle_internal_friction = MaterialModel::MaterialUtilities::phase_average_value(phase_function_values, n_phase_transitions_per_composition,
                                                                 angles_internal_friction, composition);
+            drucker_prager_parameters.angle_dilation = MaterialModel::MaterialUtilities::phase_average_value(phase_function_values, n_phase_transitions_per_composition,
+                                                       angles_dilation, composition);
             drucker_prager_parameters.cohesion = MaterialModel::MaterialUtilities::phase_average_value(phase_function_values, n_phase_transitions_per_composition,
                                                  cohesions, composition);
             drucker_prager_parameters.max_yield_stress = MaterialModel::MaterialUtilities::phase_average_value(phase_function_values, n_phase_transitions_per_composition,
@@ -220,6 +224,38 @@ namespace aspect
 
 
 
+      template <int dim>
+      std::pair<double,double>
+      DruckerPrager<dim>::compute_dilation_terms_for_stokes_system(const DruckerPragerParameters &drucker_prager_parameters,
+                                                                   const double non_yielding_viscosity,
+                                                                   const double effective_strain_rate) const
+      {
+        const double sin_phi = std::sin(drucker_prager_parameters.angle_internal_friction);
+        const double sin_psi = std::sin(drucker_prager_parameters.angle_dilation);
+
+        double alpha;
+        double alpha_bar;
+        if (dim == 2)
+          {
+            alpha     = sin_phi;
+            alpha_bar = sin_psi;
+          }
+        else if (dim == 3)
+          {
+            alpha     = 6.0 * sin_phi / (std::sqrt(3.0) * (3.0 + sin_phi));
+            alpha_bar = 6.0 * sin_psi / (std::sqrt(3.0) * (3.0 + sin_psi));
+          }
+        else
+          AssertThrow(false,ExcNotImplemented());
+
+        const double dilation_lhs_term = alpha_bar * alpha / non_yielding_viscosity;
+        const double dilation_rhs_term = alpha_bar * (2. * effective_strain_rate - drucker_prager_parameters.cohesion / non_yielding_viscosity);
+
+        return std::make_pair(dilation_lhs_term, dilation_rhs_term);
+      }
+
+
+
       template <int dim>
       void
       DruckerPrager<dim>::declare_parameters (ParameterHandler &prm)
@@ -231,6 +267,13 @@ namespace aspect
                            "those corresponding to chemical compositions. "
                            "For a value of zero, in 2d the von Mises criterion is retrieved. "
                            "Angles higher than 30 degrees are harder to solve numerically. Units: degrees.");
+        prm.declare_entry ("Angles of dilation", "0.",
+                           Patterns::Anything(),
+                           "List of angles of plastic dilation, $\\psi$, for background material and compositional fields, "
+                           "for a total of N+1 values, where N is the number of all compositional fields or only "
+                           "those corresponding to chemical compositions. "
+                           "For a value of zero, the von Mises flow rule is retrieved. "
+                           "The dilation angle should never exceed the internal friction angle.");
         prm.declare_entry ("Cohesions", "1e20",
                            Patterns::Anything(),
                            "List of cohesions, $C$, for background material and compositional fields, "
@@ -290,9 +333,26 @@ namespace aspect
         angles_internal_friction = Utilities::MapParsing::parse_map_to_double_array(prm.get("Angles of internal friction"),
                                                                                     options);
 
+        options.property_name = "Angles of dilation";
+        angles_dilation = Utilities::MapParsing::parse_map_to_double_array(prm.get("Angles of dilation"),
+                                                                           options);
+
         // Convert angles from degrees to radians
-        for (double &angle : angles_internal_friction)
-          angle *= constants::degree_to_radians;
+        for (unsigned int i = 0; i < angles_internal_friction.size(); ++i)
+          {
+            angles_internal_friction[i] *= constants::degree_to_radians;
+            angles_dilation[i]          *= constants::degree_to_radians;
+
+            AssertThrow(angles_dilation[i] <= angles_internal_friction[i],
+                        ExcMessage("The dilation angle is greater than the internal friction angle for "
+                                   "composition " + Utilities::to_string(i) + "."));
+
+            if (angles_dilation[i] > 0.0)
+              AssertThrow(this->get_parameters().enable_prescribed_dilation == true,
+                          ExcMessage("ASPECT detected a nonzero dilation angle, but dilation is "
+                                     "not enabled. Please set parameter entry 'Enable prescribed "
+                                     "dilation' to 'true'."));
+          }
 
         options.property_name = "Cohesions";
         cohesions = Utilities::MapParsing::parse_map_to_double_array(prm.get("Cohesions"),