From de158ff914dfcdd1f40bce589278b99ec7f3dd9d Mon Sep 17 00:00:00 2001 From: Juliane Dannberg Date: Sat, 15 Jun 2024 21:33:53 +0200 Subject: [PATCH] reduce grain size due to plastic strain --- include/aspect/material_model/grain_size.h | 14 ++ source/material_model/grain_size.cc | 186 +++++++++++++++------ tests/grain_size_yield.prm | 1 + tests/gs_drucker_prager.prm | 1 + 4 files changed, 153 insertions(+), 49 deletions(-) diff --git a/include/aspect/material_model/grain_size.h b/include/aspect/material_model/grain_size.h index 6aed6159ddc..63a673b3915 100644 --- a/include/aspect/material_model/grain_size.h +++ b/include/aspect/material_model/grain_size.h @@ -198,6 +198,7 @@ namespace aspect */ std::vector grain_boundary_energy; std::vector boundary_area_change_work_fraction; + std::vector yielding_work_fraction; std::vector geometric_constant; /** @@ -352,6 +353,19 @@ namespace aspect const double diffusion_viscosity, const double viscosity_guess = 0) const; + Rheology::DruckerPragerParameters compute_drucker_prager_inputs (MaterialUtilities::PhaseFunctionInputs &phase_inputs) const; + + double get_pressure_for_yielding (const double pressure, + const double adiabatic_pressure) const; + + double compute_yield_stress (const double pressure_for_yielding, + const Rheology::DruckerPragerParameters &drucker_prager_parameters) const; + + double plastic_viscosity (const double pressure_for_yielding, + const double second_strain_rate_invariant, + const double viscous_viscosity, + const Rheology::DruckerPragerParameters &drucker_prager_parameters) const; + double density (const double temperature, const double pressure, const std::vector &compositional_fields, diff --git a/source/material_model/grain_size.cc b/source/material_model/grain_size.cc index 5f807ffa1ad..28716be1c02 100644 --- a/source/material_model/grain_size.cc +++ b/source/material_model/grain_size.cc @@ -320,20 +320,47 @@ namespace aspect const double dislocation_strain_rate = second_strain_rate_invariant * current_viscosity / current_dislocation_viscosity; + double stress = 2.0 * second_strain_rate_invariant * std::min(std::max(min_eta,current_viscosity),max_eta); + + // If material is yielding, all strain rate is plastic strain rate (the rock is breaking). + bool is_yielding = false; + if (enable_drucker_prager_rheology) + { + const double gravity_norm = this->get_gravity_model().gravity_vector(in.position[i]).norm(); + const double depth = this->get_geometry_model().depth(in.position[i]); + const double rho_g = density(in.temperature[i], pressures[i], in.composition[i], in.position[i]) * gravity_norm; + + // We do not fill the phase function index, because that will be done internally in the get_phase_index() function + MaterialUtilities::PhaseFunctionInputs phase_inputs(in.temperature[i], pressures[i], depth, rho_g, numbers::invalid_unsigned_int); + + const Rheology::DruckerPragerParameters drucker_prager_parameters = compute_drucker_prager_inputs(phase_inputs); + + // We use the adiabatic pressure for grain size evolution, so we also use it for plasticity here. + const double yield_stress = compute_yield_stress(pressures[i], drucker_prager_parameters); + stress = std::min(stress, yield_stress); + is_yielding = true; + } + double grain_size_reduction_rate = 0.0; if (grain_size_evolution_formulation == Formulation::paleowattmeter) { + double strain_rate_for_reduction = boundary_area_change_work_fraction[phase_indices[i]] * dislocation_strain_rate; + if (is_yielding) + strain_rate_for_reduction = yielding_work_fraction[phase_indices[i]] * second_strain_rate_invariant; + // paleowattmeter: Austin and Evans (2007): Paleowattmeters: A scaling relation for dynamically recrystallized grain size. Geology 35, 343-346 - const double stress = 2.0 * second_strain_rate_invariant * std::min(std::max(min_eta,current_viscosity),max_eta); - grain_size_reduction_rate = 2.0 * stress * boundary_area_change_work_fraction[phase_indices[i]] * dislocation_strain_rate * grain_size * grain_size + grain_size_reduction_rate = 2.0 * stress * strain_rate_for_reduction * grain_size * grain_size / (geometric_constant[phase_indices[i]] * grain_boundary_energy[phase_indices[i]]); } else if (grain_size_evolution_formulation == Formulation::pinned_grain_damage) { + double strain_rate_for_reduction = partitioning_fraction * second_strain_rate_invariant; + if (is_yielding) + strain_rate_for_reduction = yielding_work_fraction[phase_indices[i]] * second_strain_rate_invariant; + // pinned_grain_damage: Mulyukova and Bercovici (2018) Collapse of passive margins by lithospheric damage and plunging grain size. Earth and Planetary Science Letters, 484, 341-352. - const double stress = 2.0 * second_strain_rate_invariant * std::min(std::max(min_eta,current_viscosity),max_eta); - grain_size_reduction_rate = 2.0 * stress * partitioning_fraction * second_strain_rate_invariant * grain_size * grain_size + grain_size_reduction_rate = 2.0 * stress * strain_rate_for_reduction * grain_size * grain_size * roughness_to_grain_size / (geometric_constant[phase_indices[i]] * grain_boundary_energy[phase_indices[i]] * phase_distribution); } @@ -523,6 +550,88 @@ namespace aspect + template + Rheology::DruckerPragerParameters + GrainSize:: + compute_drucker_prager_inputs (MaterialUtilities::PhaseFunctionInputs &phase_inputs) const + { + // The following handles phases + std::vector n_phases = {n_phase_transitions+1}; + std::vector phase_function_values(n_phase_transitions, 0.0); + + for (unsigned int k=0; k + double + GrainSize:: + get_pressure_for_yielding (const double pressure, + const double adiabatic_pressure) const + { + return use_adiabatic_pressure_for_yielding + ? + adiabatic_pressure + : + std::max(pressure,0.0); + } + + + template + double + GrainSize:: + compute_yield_stress (const double pressure_for_yielding, + const Rheology::DruckerPragerParameters &drucker_prager_parameters) const + { + return drucker_prager_plasticity.compute_yield_stress(drucker_prager_parameters.cohesion, + drucker_prager_parameters.angle_internal_friction, + pressure_for_yielding, + drucker_prager_parameters.max_yield_stress); + } + + + + template + double + GrainSize:: + plastic_viscosity (const double pressure_for_yielding, + const double second_strain_rate_invariant, + const double viscous_viscosity, + const Rheology::DruckerPragerParameters &drucker_prager_parameters) const + { + // Calculate non-yielding (viscous) stress magnitude. + const double non_yielding_stress = 2. * viscous_viscosity * second_strain_rate_invariant; + const double yield_stress = compute_yield_stress (pressure_for_yielding, drucker_prager_parameters); + + double effective_viscosity = viscous_viscosity; + + // Apply plastic yielding: + // If the non-yielding stress is greater than the yield stress, + // rescale the viscosity back to yield surface + if (non_yielding_stress >= yield_stress) + { + effective_viscosity = drucker_prager_plasticity.compute_viscosity(drucker_prager_parameters.cohesion, + drucker_prager_parameters.angle_internal_friction, + pressure_for_yielding, + second_strain_rate_invariant, + drucker_prager_parameters.max_yield_stress, + effective_viscosity); + } + return effective_viscosity; + } + + + template double GrainSize:: @@ -817,8 +926,6 @@ namespace aspect if (in.requests_property(MaterialProperties::viscosity)) { - double effective_viscosity; - double disl_viscosity = std::numeric_limits::max(); Assert(std::isfinite(in.strain_rate[i].norm()), ExcMessage("Invalid strain_rate in the MaterialModelInputs. This is likely because it was " "not filled by the caller.")); @@ -842,62 +949,31 @@ namespace aspect second_strain_rate_invariant, phase_indices[i]); + double effective_viscosity = diff_viscosity; + double disl_viscosity = std::numeric_limits::max(); + if (std::abs(second_strain_rate_invariant) > 1e-30) { disl_viscosity = dislocation_viscosity(in.temperature[i], adiabatic_temperature, adiabatic_pressures[i], in.strain_rate[i], phase_indices[i], diff_viscosity); effective_viscosity = disl_viscosity * diff_viscosity / (disl_viscosity + diff_viscosity); } - else - effective_viscosity = diff_viscosity; + const double non_yielding_stress = 2. * effective_viscosity * second_strain_rate_invariant; if (enable_drucker_prager_rheology) { - // Calculate non-yielding (viscous) stress magnitude. - const double non_yielding_stress = 2. * effective_viscosity * second_strain_rate_invariant; - - // The following handles phases - std::vector n_phases = {n_phase_transitions+1}; - std::vector phase_function_values(n_phase_transitions, 0.0); - - for (unsigned int k=0; k= yield_stress) - { - effective_viscosity = drucker_prager_plasticity.compute_viscosity(drucker_prager_parameters.cohesion, - drucker_prager_parameters.angle_internal_friction, - pressure_for_yielding, - second_strain_rate_invariant, - drucker_prager_parameters.max_yield_stress, - effective_viscosity); - } + const Rheology::DruckerPragerParameters drucker_prager_parameters = compute_drucker_prager_inputs(phase_inputs); + const double pressure_for_yielding = get_pressure_for_yielding (in.pressure[i], adiabatic_pressures[i]); + effective_viscosity = plastic_viscosity (pressure_for_yielding, + second_strain_rate_invariant, + effective_viscosity, + drucker_prager_parameters); PlasticAdditionalOutputs *plastic_out = out.template get_additional_output>(); if (plastic_out != nullptr) { + const double yield_stress = compute_yield_stress(pressure_for_yielding, drucker_prager_parameters); + plastic_out->cohesions[i] = drucker_prager_parameters.cohesion; plastic_out->friction_angles[i] = drucker_prager_parameters.angle_internal_friction; plastic_out->yield_stresses[i] = yield_stress; @@ -1113,6 +1189,11 @@ namespace aspect "The fraction $\\chi$ of work done by dislocation creep to change the grain boundary area. " "List must have one more entry than the Phase transition depths. " "Units: \\si{\\joule\\per\\meter\\squared}."); + prm.declare_entry ("Work fraction for boundary area change when yielding", "0.2", + Patterns::List (Patterns::Double (0.)), + "The fraction $\\chi$ of work done by plastic tielding to change the grain boundary area. " + "List must have one more entry than the Phase transition depths. " + "Units: \\si{\\joule\\per\\meter\\squared}."); prm.declare_entry ("Geometric constant", "3.", Patterns::List (Patterns::Double (0.)), "The geometric constant $c$ used in the paleowattmeter grain size reduction law. " @@ -1396,6 +1477,8 @@ namespace aspect (Utilities::split_string_list(prm.get ("Average specific grain boundary energy"))); boundary_area_change_work_fraction = Utilities::string_to_double (Utilities::split_string_list(prm.get ("Work fraction for boundary area change"))); + yielding_work_fraction = Utilities::string_to_double + (Utilities::split_string_list(prm.get ("Work fraction for boundary area change when yielding"))); geometric_constant = Utilities::string_to_double (Utilities::split_string_list(prm.get ("Geometric constant"))); @@ -1567,6 +1650,11 @@ namespace aspect std::vector n_phases = {n_phase_transitions+1}; drucker_prager_plasticity.parse_parameters(prm, std::make_unique> (n_phases)); + + if (enable_drucker_prager_rheology) + AssertThrow(n_phases[0] == yielding_work_fraction.size(), + ExcMessage("Error: The list of ``Work fraction for boundary area change when yielding'' " + "does not have the correct length! It needs to have one entry per phase.")); } prm.leave_subsection(); } diff --git a/tests/grain_size_yield.prm b/tests/grain_size_yield.prm index 4a8e5c04db4..361b2c6d44f 100644 --- a/tests/grain_size_yield.prm +++ b/tests/grain_size_yield.prm @@ -19,6 +19,7 @@ subsection Material model set Reference temperature = 293 set Grain growth rate constant = 0 set Work fraction for boundary area change = 0 + set Work fraction for boundary area change when yielding = 0 # Diffusion creep # new scaled prefactors to match vertical viscosity profile diff --git a/tests/gs_drucker_prager.prm b/tests/gs_drucker_prager.prm index 311c133c4d1..db07d027fcf 100644 --- a/tests/gs_drucker_prager.prm +++ b/tests/gs_drucker_prager.prm @@ -90,6 +90,7 @@ subsection Material model set Reference temperature = 1600 set Grain growth rate constant = 0 set Work fraction for boundary area change = 0 + set Work fraction for boundary area change when yielding = 0 # Faul and Jackson 2007 # Diffusion creep