Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Newton method for plastic dilation #5953

Open
wants to merge 4 commits into
base: main
Choose a base branch
from
Open

Newton method for plastic dilation #5953

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
6463b7f
Add dilation relevant terms to the Newton solver
Jun 30, 2024
b0e8fa4
fix bugs: missing pressure_scaling in local_apply()
Jul 25, 2024
251f144
change the eos back to incompressible model
Jul 25, 2024
b1c748e
treat plastic dilation in the same way as compressible model
Jul 25, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
17 changes: 17 additions & 0 deletions include/aspect/material_model/rheology/drucker_prager.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,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
*/
Expand Down Expand Up @@ -138,9 +143,21 @@ namespace aspect
compute_derivative (const double angle_internal_friction,
const double effective_strain_rate) const;

/**
* Compute the plastic dilation term, which is defined as the product of
* the plastic multiplier and the derivative of plastic potential with
* respect to minus pressure.
*/
double
compute_plastic_dilation (const double angle_dilation,
const double non_yielding_viscosity,
const double effective_viscosity,
const double effective_strain_rate) const;

private:

std::vector<double> angles_internal_friction;
std::vector<double> angles_dilation;
std::vector<double> cohesions;
double max_yield_stress;

Expand Down
9 changes: 8 additions & 1 deletion include/aspect/material_model/rheology/visco_plastic.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -108,6 +108,13 @@ namespace aspect
* The current cohesion.
*/
std::vector<double> current_cohesions;

/**
* The plastic dilation terms, defined as the product of the
* plastic multiplier and the derivative of plastic potential
* with respect to negative pressure.
*/
std::vector<double> plastic_dilation;
};

namespace Rheology
Expand Down Expand Up @@ -149,7 +156,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>(),
Expand Down
10 changes: 10 additions & 0 deletions include/aspect/newton.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,6 +63,16 @@ namespace aspect
* derivatives when material averaging is applied.
*/
std::vector<double> viscosity_derivative_averaging_weights;

/**
* The derivatives of the prescribed dilation with respect to strain rate.
*/
std::vector<SymmetricTensor<2,dim>> dilation_derivative_wrt_strain_rate;

/**
* The derivatives of the prescribed dilation term with respect to pressure.
*/
std::vector<double> dilation_derivative_wrt_pressure;
};
}

Expand Down
23 changes: 21 additions & 2 deletions include/aspect/stokes_matrix_free.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -125,6 +125,12 @@ namespace aspect
*/
bool apply_stabilization_free_surface_faces;

/**
* If true, derivatives of the prescribed dilation are
* part of the operator.
*/
bool enable_prescribed_dilation;

/**
* Table which stores viscosity values for each cell.
*
Expand All @@ -144,7 +150,7 @@ namespace aspect
* variables: viscosity derivative with respect to pressure,
* the Newton derivative scaling factor, and the averaging weight.
*/
Table<2, VectorizedArray<number>> newton_factor_wrt_pressure_table;
Table<2, VectorizedArray<number>> viscosity_derivative_wrt_pressure_table;

/**
* Table which stores the product of the following four
Expand All @@ -154,7 +160,20 @@ namespace aspect
* is PD or SPD, otherwise, it is 1.
*/
Table<2, SymmetricTensor<2, dim, VectorizedArray<number>>>
newton_factor_wrt_strain_rate_table;
viscosity_derivative_wrt_strain_rate_table;

/**
* Table which stores the product of the dilation derivative
* with respect to pressure and the Newton derivative scaling factor.
*/
Table<2, VectorizedArray<number>> dilation_derivative_wrt_pressure_table;

/**
* Table which stores the product of the 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
Expand Down
50 changes: 48 additions & 2 deletions source/material_model/rheology/drucker_prager.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,7 @@ namespace aspect
{
DruckerPragerParameters::DruckerPragerParameters()
: angle_internal_friction (numbers::signaling_nan<double>()),
angle_dilation (numbers::signaling_nan<double>()),
cohesion (numbers::signaling_nan<double>()),
max_yield_stress (numbers::signaling_nan<double>())
{}
Expand All @@ -60,13 +61,16 @@ 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];
}
else
{
// 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);
}
Expand Down Expand Up @@ -181,6 +185,26 @@ namespace aspect



template <int dim>
double
DruckerPrager<dim>::compute_plastic_dilation(const double angle_dilation,
const double non_yielding_viscosity,
const double effective_viscosity,
const double effective_strain_rate) const
{
const double sin_psi = std::sin(angle_dilation);
const double minus_dQ_dp = (dim == 3 ?
6.0 * sin_psi / (std::sqrt(3.0) * (3.0 + sin_psi)) :
sin_psi);

const double plastic_multiplier =
std::max(0.0, (1.0 - effective_viscosity / non_yielding_viscosity) * 2.0 * effective_strain_rate);

return minus_dQ_dp * plastic_multiplier;
}



template <int dim>
void
DruckerPrager<dim>::declare_parameters (ParameterHandler &prm)
Expand All @@ -192,6 +216,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, "
Expand Down Expand Up @@ -248,10 +279,25 @@ namespace aspect

angles_internal_friction = Utilities::MapParsing::parse_map_to_double_array(prm.get("Angles of internal friction"),
options);
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"),
Expand Down
Loading