From b1e55a069c0234abee7673aa7515e0a254a99689 Mon Sep 17 00:00:00 2001 From: Bob Myhill Date: Thu, 25 Jul 2024 16:58:15 +0100 Subject: [PATCH] added isostress averaging --- .../rheology/composite_visco_plastic.h | 37 ++ .../rheology/composite_visco_plastic.cc | 315 ++++++++++++++++-- 2 files changed, 333 insertions(+), 19 deletions(-) diff --git a/include/aspect/material_model/rheology/composite_visco_plastic.h b/include/aspect/material_model/rheology/composite_visco_plastic.h index 02b6e107b92..6ee3a6c764d 100644 --- a/include/aspect/material_model/rheology/composite_visco_plastic.h +++ b/include/aspect/material_model/rheology/composite_visco_plastic.h @@ -83,6 +83,25 @@ namespace aspect const std::vector &phase_function_values = std::vector(), const std::vector &n_phase_transitions_per_composition = std::vector()) const; + /** + * Compute the isostress viscosity over all compositional fields + * based on the composite viscous creep law. + * If @p expected_n_phases_per_composition points to a vector of + * unsigned integers this is considered the number of phase transitions + * for each compositional field and viscosity will be first computed on + * each phase and then averaged for each compositional field. + */ + double + compute_isostress_viscosity (const double pressure, + const double temperature, + const double grain_size, + const std::vector &volume_fractions, + const SymmetricTensor<2,dim> &strain_rate, + std::vector &partial_strain_rates, + const std::vector &phase_function_values = std::vector(), + const std::vector &n_phase_transitions_per_composition = std::vector()) const; + + /** * Compute the compositional field viscosity * based on the composite viscous creep law. @@ -111,7 +130,25 @@ namespace aspect const double viscoplastic_stress, std::vector &partial_strain_rates) const; + /** + * Compute the total strain rate and the first derivative of log strain rate + * with respect to log viscoplastic stress from individual log strain rate components. + * Also updates the partial_strain_rates vector + */ + std::pair + calculate_log_strain_rate_and_derivative(const std::vector, 4>> &logarithmic_strain_rates_and_stress_derivatives, + const double viscoplastic_stress, + std::vector &partial_strain_rates) const; + private: + /** + * Enumeration for selecting which type of viscosity averaging to use. + */ + enum ViscosityAveragingScheme + { + isostrain, + isostress + } viscosity_averaging_scheme; /** * Whether to use different deformation mechanisms diff --git a/source/material_model/rheology/composite_visco_plastic.cc b/source/material_model/rheology/composite_visco_plastic.cc index d8537ca5aa8..5ee4f4e2acf 100644 --- a/source/material_model/rheology/composite_visco_plastic.cc +++ b/source/material_model/rheology/composite_visco_plastic.cc @@ -82,38 +82,243 @@ namespace aspect double viscosity = 0.; partial_strain_rates.resize(5, 0.); + // Isostress averaging + if (viscosity_averaging_scheme == isostress) + { + viscosity += compute_isostress_viscosity (pressure, + temperature, + grain_size, + volume_fractions, + strain_rate, + partial_strain_rates, + phase_function_values, + n_phase_transitions_per_composition); + } + // Isostrain averaging + if (viscosity_averaging_scheme == isostrain) + { + double total_volume_fraction = 1.; + for (unsigned int composition=0; composition < number_of_chemical_compositions; ++composition) + { + // Only include the contribution to the viscosity + // from a given composition if the volume fraction exceeds + // a certain (small) fraction. + if (volume_fractions[composition] > 2.*std::numeric_limits::epsilon()) + { + std::vector partial_strain_rates_composition(5, 0.); + viscosity += (volume_fractions[composition] + * compute_composition_viscosity (pressure, + temperature, + grain_size, + composition, + strain_rate, + partial_strain_rates_composition, + phase_function_values, + n_phase_transitions_per_composition)); + for (unsigned int j=0; j < 5; ++j) + partial_strain_rates[j] += volume_fractions[composition] * partial_strain_rates_composition[j]; + } + else + { + total_volume_fraction -= volume_fractions[composition]; + } + + viscosity /= total_volume_fraction; + for (unsigned int j=0; j < 5; ++j) + partial_strain_rates[j] /= total_volume_fraction; + } + } + return viscosity; + } + + + + + template + double + CompositeViscoPlastic::compute_isostress_viscosity (const double pressure, + const double temperature, + const double grain_size, + const std::vector &volume_fractions, + const SymmetricTensor<2,dim> &strain_rate, + std::vector &partial_strain_rates, + const std::vector &phase_function_values, + const std::vector &n_phase_transitions_per_composition) const + { + // If strain rate is zero (like during the first time step) set it to some very small number + // to prevent a division-by-zero, and a floating point exception. + // Otherwise, calculate the square-root of the norm of the second invariant of the deviatoric- + // strain rate (often simplified as epsilondot_ii) + const double edot_ii = std::max(std::sqrt(std::max(-second_invariant(deviator(strain_rate)), 0.)), + minimum_strain_rate); + const double log_edot_ii = std::log(edot_ii); + + std::vector diffusion_creep_parameters; + std::vector dislocation_creep_parameters; + std::vector peierls_creep_parameters; + std::vector drucker_prager_parameters; + + // 1) Estimate the stress running through the creep elements and the + // maximum viscosity element, whose viscosity is defined in parse_parameters. + // The creep elements for each composition are arranged in series. + // Taking the minimum viscosity from + // all of these elements provides an excellent first approximation + // to the true viscosity of that composition. + // For a starting guess, we further assume that one composition takes + // up all the strain. + + // The stress can then be calculated as 2 * viscoplastic_viscosity_guess * edot_ii + double viscoplastic_viscosity_guess = maximum_viscosity; + double total_volume_fraction = 1.; - for (unsigned int composition=0; composition < number_of_chemical_compositions; ++composition) + std::vector active_compositions; + for (unsigned int composition = 0; composition < number_of_chemical_compositions; ++composition) { // Only include the contribution to the viscosity // from a given composition if the volume fraction exceeds // a certain (small) fraction. - if (volume_fractions[composition] > 2.*std::numeric_limits::epsilon()) + const double volume_fraction = volume_fractions[composition]; + if (volume_fraction > 2. * std::numeric_limits::epsilon()) { - std::vector partial_strain_rates_composition(5, 0.); - viscosity += (volume_fractions[composition] - * compute_composition_viscosity (pressure, - temperature, - grain_size, - composition, - strain_rate, - partial_strain_rates_composition, - phase_function_values, - n_phase_transitions_per_composition)); - for (unsigned int j=0; j < 5; ++j) - partial_strain_rates[j] += volume_fractions[composition] * partial_strain_rates_composition[j]; + active_compositions.push_back(composition); + + if (use_diffusion_creep) + { + diffusion_creep_parameters.push_back(diffusion_creep->compute_creep_parameters(composition, phase_function_values, n_phase_transitions_per_composition)); + viscoplastic_viscosity_guess = std::min(diffusion_creep->compute_viscosity(pressure, temperature, grain_size, composition, phase_function_values, n_phase_transitions_per_composition)/volume_fraction, viscoplastic_viscosity_guess); + } + + if (use_dislocation_creep) + { + dislocation_creep_parameters.push_back(dislocation_creep->compute_creep_parameters(composition, phase_function_values, n_phase_transitions_per_composition)); + viscoplastic_viscosity_guess = std::min(dislocation_creep->compute_viscosity(edot_ii/volume_fraction, pressure, temperature, composition, phase_function_values, n_phase_transitions_per_composition)/volume_fraction, viscoplastic_viscosity_guess); + } + + if (use_peierls_creep) + { + peierls_creep_parameters.push_back(peierls_creep->compute_creep_parameters(composition)); + viscoplastic_viscosity_guess = std::min(peierls_creep->compute_approximate_viscosity(edot_ii/volume_fraction, pressure, temperature, composition)/volume_fraction, viscoplastic_viscosity_guess); + } + + if (use_drucker_prager) + { + Rheology::DruckerPragerParameters drucker_prager_parameters_c = drucker_prager->compute_drucker_prager_parameters(composition, phase_function_values, n_phase_transitions_per_composition); + drucker_prager_parameters.push_back(drucker_prager_parameters_c); + viscoplastic_viscosity_guess = std::min(drucker_prager->compute_viscosity(drucker_prager_parameters_c.cohesion, + drucker_prager_parameters_c.angle_internal_friction, pressure, edot_ii/volume_fraction, drucker_prager_parameters_c.max_yield_stress)/volume_fraction, + viscoplastic_viscosity_guess); + } } else { total_volume_fraction -= volume_fractions[composition]; } + } + + double viscoplastic_stress = 2. * viscoplastic_viscosity_guess * edot_ii / total_volume_fraction; + double log_viscoplastic_stress = std::log(viscoplastic_stress); + + // 2) Calculate the log strain rate and first derivative with respect to + // log stress for each element using the guessed creep stress. + std::vector, 4>> log_edot_and_deriv(active_compositions.size()); - viscosity /= total_volume_fraction; - for (unsigned int j=0; j < 5; ++j) - partial_strain_rates[j] /= total_volume_fraction; + for (unsigned int i = 0; i < active_compositions.size(); ++i) + { + const double log_volume_fraction = std::log(volume_fractions[active_compositions[i]]); + + if (use_diffusion_creep) + { + log_edot_and_deriv[i][0] = diffusion_creep->compute_log_strain_rate_and_derivative(log_viscoplastic_stress, pressure, temperature, grain_size, diffusion_creep_parameters[i]); + log_edot_and_deriv[i][0].first += log_volume_fraction; + } + if (use_dislocation_creep) + { + log_edot_and_deriv[i][1] = dislocation_creep->compute_log_strain_rate_and_derivative(log_viscoplastic_stress, pressure, temperature, dislocation_creep_parameters[i]); + log_edot_and_deriv[i][1].first += log_volume_fraction; + } + if (use_peierls_creep) + { + log_edot_and_deriv[i][2] = peierls_creep->compute_approximate_log_strain_rate_and_derivative(log_viscoplastic_stress, pressure, temperature, peierls_creep_parameters[i]); + log_edot_and_deriv[i][2].first += log_volume_fraction; + } + if (use_drucker_prager) + { + log_edot_and_deriv[i][3] = drucker_prager->compute_log_strain_rate_and_derivative(log_viscoplastic_stress, pressure, drucker_prager_parameters[i]); + log_edot_and_deriv[i][3].first += log_volume_fraction; + } } - return viscosity; + + // 3) Calculate the total log strain rate from the first estimates + // for the component log strain rates. + // This will generally be more than the input total strain rate because + // the creep stress in Step 1 was been calculated assuming that + // only one mechanism was active, whereas the strain rate + // calculated in Step 2 allowed all the mechanisms to + // accommodate strain at that creep stress. + std::pair log_edot_ii_and_deriv_iterate = calculate_log_strain_rate_and_derivative(log_edot_and_deriv, + viscoplastic_stress, + partial_strain_rates); + double log_strain_rate_residual = log_edot_ii_and_deriv_iterate.first - log_edot_ii; + + // 4) In this rheology model, the total strain rate is partitioned between + // different flow components. We do not know how the strain is partitioned + // between these components. + + // The following while loop contains a Newton iteration to obtain the + // viscoplastic stress that is consistent with the total strain rate. + unsigned int stress_iteration = 0; + while (std::abs(log_strain_rate_residual) > log_strain_rate_residual_threshold && stress_iteration < stress_max_iteration_number) + { + // Apply the Newton update for the log creep stress using the + // strain-rate residual and strain-rate stress derivative + double delta_log_viscoplastic_stress = log_strain_rate_residual / log_edot_ii_and_deriv_iterate.second; + log_viscoplastic_stress += delta_log_viscoplastic_stress; + viscoplastic_stress = std::exp(log_viscoplastic_stress); + + // Update the strain rates of all mechanisms with the new stress + for (auto &log_edot_and_deriv_c : log_edot_and_deriv) + { + for (auto &i : active_flow_mechanisms) + log_edot_and_deriv_c[i].first += log_edot_and_deriv_c[i].second * delta_log_viscoplastic_stress; + } + + // Compute the new log strain rate residual and log stress derivative + log_edot_ii_and_deriv_iterate = calculate_log_strain_rate_and_derivative(log_edot_and_deriv, + viscoplastic_stress, + partial_strain_rates); + log_strain_rate_residual = log_edot_ii - log_edot_ii_and_deriv_iterate.first; + + ++stress_iteration; + + // If anything that would be used in the next iteration is not finite, the + // Newton iteration would trigger an exception and we want to abort the + // iteration instead. + // Currently, we still throw an exception, but if this exception is thrown, + // another more robust iterative scheme should be implemented + // (similar to that seen in the diffusion-dislocation material model). + const bool abort_newton_iteration = !numbers::is_finite(log_viscoplastic_stress) || !numbers::is_finite(log_strain_rate_residual) || stress_iteration == stress_max_iteration_number; + AssertThrow(!abort_newton_iteration, + ExcMessage("No convergence has been reached in the loop that determines " + "the composite viscous creep stress. Aborting! " + "Residual is " + + Utilities::to_string(log_strain_rate_residual) + + " after " + Utilities::to_string(stress_iteration) + " iterations. " + "You can increase the number of iterations by adapting the " + "parameter 'Maximum creep strain rate iterations'.")); + } + + // 5) We have now calculated the viscoplastic stress consistent with the total + // strain rate, but the viscoplastic stress is only one component of the total stress, + // because this material model also includes a viscosity damper + // arranged in parallel with the viscoplastic elements. + // The total stress is equal to the sum of the viscoplastic stress and + // minimum stress. + const double damper_stress = 2. * damper_viscosity * edot_ii; + const double total_stress = viscoplastic_stress + damper_stress; + + // 6) Return the effective creep viscosity using the total stress + return total_stress / (2. * edot_ii); } @@ -311,7 +516,63 @@ namespace aspect const double log_viscoplastic_strain_rate_derivative = weighted_stress_derivative_sum / viscoplastic_strain_rate_sum; - // Some opaque mathmatics converts the viscoplastic strain rate to the total strain rate. + // Some opaque mathematics converts the viscoplastic strain rate to the total strain rate. + const double f = viscoplastic_stress / (2. * maximum_viscosity); + const double strain_rate = (strain_rate_scaling_factor * viscoplastic_strain_rate_sum) + f; + partial_strain_rates[4] = strain_rate - viscoplastic_strain_rate_sum; + // And the partial derivative of the log *total* strain rate + // with respect to log *viscoplastic* stress follows as + const double log_strain_rate_derivative = (strain_rate_scaling_factor * viscoplastic_strain_rate_sum * log_viscoplastic_strain_rate_derivative + f) / strain_rate; + + return std::make_pair(std::log(strain_rate), log_strain_rate_derivative); + } + + + + template + std::pair + CompositeViscoPlastic::calculate_log_strain_rate_and_derivative(const std::vector, 4>> &logarithmic_strain_rates_and_stress_derivatives, + const double viscoplastic_stress, + std::vector &partial_strain_rates) const + { + // The total strain rate + double viscoplastic_strain_rate_sum = 0.0; + + // The sum of the stress derivatives multiplied by the mechanism strain rates + double weighted_stress_derivative_sum = 0.0; + + // The first derivative of log(strain rate) with respect to log(stress) + // is computed as sum_i(stress_exponent_i * edot_i) / sum_i(edot_i) + // i.e., the stress exponents weighted by strain rate fraction + // summed over the individual flow mechanisms (i). + // Loop over active flow laws and add their contributions + // to the strain rate and stress derivative + + // First, make sure that the active partial_strain_rates are reset to zero. + for (auto &i : active_flow_mechanisms) + partial_strain_rates[i] = 0; + + for (auto &logarithmic_strain_rates_and_stress_derivatives_c : logarithmic_strain_rates_and_stress_derivatives) + { + for (auto &i : active_flow_mechanisms) + { + double mechanism_log_strain_rate = logarithmic_strain_rates_and_stress_derivatives_c[i].first; + + // Check if the mechanism strain rate is within bounds to prevent underflow + if (mechanism_log_strain_rate >= logmin) + { + const double mechanism_strain_rate = std::exp(mechanism_log_strain_rate); + partial_strain_rates[i] += mechanism_strain_rate; + const double log_stress_derivative = logarithmic_strain_rates_and_stress_derivatives_c[i].second; + viscoplastic_strain_rate_sum += mechanism_strain_rate; + weighted_stress_derivative_sum += log_stress_derivative * mechanism_strain_rate; + } + } + } + + const double log_viscoplastic_strain_rate_derivative = weighted_stress_derivative_sum / viscoplastic_strain_rate_sum; + + // Some opaque mathematics converts the viscoplastic strain rate to the total strain rate. const double f = viscoplastic_stress / (2. * maximum_viscosity); const double strain_rate = (strain_rate_scaling_factor * viscoplastic_strain_rate_sum) + f; partial_strain_rates[4] = strain_rate - viscoplastic_strain_rate_sum; @@ -328,6 +589,13 @@ namespace aspect void CompositeViscoPlastic::declare_parameters (ParameterHandler &prm) { + prm.declare_entry ("Viscosity averaging scheme", "isostress", + Patterns::Selection("isostress|isostrain|unspecified"), + "Determines the relationship between the conditions experienced by " + "each chemical compositional field in a composite material. " + "Select either isostress (the default option, Reuss averaging) " + "or isostrain (Voigt averaging)."); + prm.declare_entry ("Include diffusion creep in composite rheology", "true", Patterns::Bool (), "Whether to include diffusion creep in the composite rheology formulation."); @@ -390,6 +658,15 @@ namespace aspect CompositeViscoPlastic::parse_parameters (ParameterHandler &prm, const std::unique_ptr> &expected_n_phases_per_composition) { + if (prm.get ("Viscosity averaging scheme") == "isostress") + viscosity_averaging_scheme = isostress; + else if (prm.get ("Viscosity averaging scheme") == "isostrain") + viscosity_averaging_scheme = isostrain; + else + { + AssertThrow(false, ExcMessage("Not a valid viscosity averaging scheme. Choose either isostress or isostrain.")); + } + // A background field is required by the subordinate material models number_of_chemical_compositions = this->introspection().n_chemical_composition_fields() + 1;