From 561970e4a16314589717c97767b074965502b22f Mon Sep 17 00:00:00 2001 From: Bob Myhill Date: Wed, 9 Dec 2020 16:10:21 +0000 Subject: [PATCH] add volume fraction limit for viscosity calculation --- .../rheology/composite_visco_plastic.cc | 91 ++++++++++++------- 1 file changed, 58 insertions(+), 33 deletions(-) diff --git a/source/material_model/rheology/composite_visco_plastic.cc b/source/material_model/rheology/composite_visco_plastic.cc index 68a4b8e9d12..340e1107475 100644 --- a/source/material_model/rheology/composite_visco_plastic.cc +++ b/source/material_model/rheology/composite_visco_plastic.cc @@ -120,11 +120,10 @@ namespace aspect { 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; } + viscosity /= total_volume_fraction; + for (unsigned int j=0; j < 5; ++j) + partial_strain_rates[j] /= total_volume_fraction; } // Final step for Voigt-Reuss-Hill averaging is to divide @@ -220,22 +219,35 @@ namespace aspect while (std::abs(strain_rate_residual) > strain_rate_residual_threshold && stress_iteration < stress_max_iteration_number) { - + double total_volume_fraction = 1.; std::pair creep_edot_and_deriv = std::make_pair(0., 0.); for (unsigned int composition=0; composition < number_of_compositions; ++composition) { - creep_edot_and_deriv = (creep_edot_and_deriv - + volume_fractions[composition] - * compute_strain_rate_and_derivative (creep_stress, - pressure, - temperature, - composition, - diffusion_creep_parameters[composition], - dislocation_creep_parameters[composition], - peierls_creep_parameters[composition], - drucker_prager_parameters)); + // 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()) + { + creep_edot_and_deriv = (creep_edot_and_deriv + + volume_fractions[composition] + * compute_strain_rate_and_derivative (creep_stress, + pressure, + temperature, + composition, + diffusion_creep_parameters[composition], + dislocation_creep_parameters[composition], + peierls_creep_parameters[composition], + drucker_prager_parameters)); + } + else + { + total_volume_fraction -= volume_fractions[composition]; + } } + creep_edot_and_deriv = std::make_pair(creep_edot_and_deriv.first/total_volume_fraction, + creep_edot_and_deriv.second/total_volume_fraction); + const double strain_rate = creep_stress/(2.*max_viscosity) + (max_viscosity/(max_viscosity - min_viscosity))*creep_edot_and_deriv.first; strain_rate_deriv = 1./(2.*max_viscosity) + (max_viscosity/(max_viscosity - min_viscosity))*creep_edot_and_deriv.second; @@ -277,35 +289,48 @@ namespace aspect // The components of partial_strain_rates must be provided in the order // dictated by make_strain_rate_additional_outputs_names + double total_volume_fraction = 1.; std::fill(partial_strain_rates.begin(), partial_strain_rates.end(), 0); for (unsigned int composition=0; composition < number_of_compositions; ++composition) { - if (use_diffusion_creep) + // 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 std::pair diff_edot_and_deriv = diffusion_creep->compute_strain_rate_and_derivative(creep_stress, pressure, temperature, diffusion_creep_parameters[composition]); - partial_strain_rates[0] = volume_fractions[composition] * diff_edot_and_deriv.first; - } + if (use_diffusion_creep) + { + const std::pair diff_edot_and_deriv = diffusion_creep->compute_strain_rate_and_derivative(creep_stress, pressure, temperature, diffusion_creep_parameters[composition]); + partial_strain_rates[0] = volume_fractions[composition] * diff_edot_and_deriv.first; + } - if (use_dislocation_creep) - { - const std::pair disl_edot_and_deriv = dislocation_creep->compute_strain_rate_and_derivative(creep_stress, pressure, temperature, dislocation_creep_parameters[composition]); - partial_strain_rates[1] = volume_fractions[composition] * disl_edot_and_deriv.first; - } + if (use_dislocation_creep) + { + const std::pair disl_edot_and_deriv = dislocation_creep->compute_strain_rate_and_derivative(creep_stress, pressure, temperature, dislocation_creep_parameters[composition]); + partial_strain_rates[1] = volume_fractions[composition] * disl_edot_and_deriv.first; + } - if (use_peierls_creep) - { - const std::pair prls_edot_and_deriv = peierls_creep->compute_strain_rate_and_derivative(creep_stress, pressure, temperature, peierls_creep_parameters[composition]); - partial_strain_rates[2] = volume_fractions[composition] * prls_edot_and_deriv.first; - } + if (use_peierls_creep) + { + const std::pair prls_edot_and_deriv = peierls_creep->compute_strain_rate_and_derivative(creep_stress, pressure, temperature, peierls_creep_parameters[composition]); + partial_strain_rates[2] = volume_fractions[composition] * prls_edot_and_deriv.first; + } - if (use_drucker_prager) + if (use_drucker_prager) + { + const std::pair drpr_edot_and_deriv = drucker_prager->compute_strain_rate_and_derivative(creep_stress, pressure, composition, drucker_prager_parameters); + partial_strain_rates[3] = volume_fractions[composition] * drpr_edot_and_deriv.first; + } + } + else { - const std::pair drpr_edot_and_deriv = drucker_prager->compute_strain_rate_and_derivative(creep_stress, pressure, composition, drucker_prager_parameters); - partial_strain_rates[3] = volume_fractions[composition] * drpr_edot_and_deriv.first; + total_volume_fraction -= volume_fractions[composition]; } - } + for (unsigned int j=0; j < 4; ++j) + partial_strain_rates[j] /= total_volume_fraction; + partial_strain_rates[4] = total_stress/(2.*max_viscosity); // Now we return the viscosity using the total stress