Skip to content

Add temperature parameters for RefJC viscoplasticity, along with GP output #1879

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

Open
dragos-ana wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Add temperature parameters for RefJC viscoplasticity, along with GP output #1879

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
36 changes: 29 additions & 7 deletions src/global_legacy_module/4C_global_legacy_module_validmaterials.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2873,23 +2873,45 @@ std::unordered_map<Core::Materials::MaterialType, Core::IO::InputSpec> Global::v

/*----------------------------------------------------------------------*/
{
using namespace Core::IO::InputSpecBuilders::Validators;

known_materials[Core::Materials::mvl_reformulated_Johnson_Cook] =
group("MAT_ViscoplasticLawReformulatedJohnsonCook",
{
parameter<double>("STRAIN_RATE_PREFAC",
{.description = "reference plastic strain rate $\\dot{P}_0$ "}),
{.description = "reference plastic strain rate $\\dot{P}_0$",
.validator = positive<double>()}),
parameter<double>("STRAIN_RATE_EXP_FAC",
{.description = "exponential factor of plastic strain rate $C$"}),
{.description = "exponential factor of plastic strain rate $C$",
.validator = positive<double>()}),
parameter<double>("INIT_YIELD_STRENGTH",
{.description = "initial yield strength of the material $A_0$"}),
{.description = "initial yield strength of the material $A_0$",
.validator = positive<double>()}),
parameter<double>("ISOTROP_HARDEN_PREFAC",
{.description = "prefactor of the isotropic hardening stress $B_0$"}),
{.description =
"prefactor of the isotropic hardening stress / hardening modulus $B_0$",
.validator = positive_or_zero<double>()}),
parameter<double>("ISOTROP_HARDEN_EXP",
{.description = "exponent of the isotropic hardening stress $n$"}),
{.description = "exponent of the isotropic hardening stress $n$",
.validator = positive_or_zero<double>()}),
parameter<double>("REF_TEMPERATURE",
{.description = "reference temperature $T_0$ for evaluating the "
"yield strength $A_0$ and the hardening "
"modulus $B_0$",
.validator = positive<double>()}),
parameter<double>(
"MELT_TEMPERATURE", {.description = "melting temperature $T_{\\mathrm{melt}}$",
.validator = positive<double>()}),
parameter<double>("TEMPERATURE_SENS", {.description = "temperature sensitivity $m$",
.validator = positive_or_zero<double>()}),

},
{.description = "Reformulation of the Johnson-Cook viscoplastic law (comprising flow "
"rule $\\dot{P} = \\dot{P}_0 \\exp \\left( \\frac{ \\Sigma_{eq}}{C "
"\\Sigma_y} - \\frac{1}{C} \\right) - \\dot{P}_0$ and hardening law), "
"rule $\\dot{P} = \\dot{P}_0 \\exp \\left( \\frac{ \\sigma_{eq}}{C "
"\\sigma_{\\mathrm{Y}}} - \\frac{1}{C} \\right) - \\dot{P}_0$ and "
"hardening law $\\sigma_{\\mathrm{Y}} = (A_0 + "
"B_0 \\cdot P^{n}) \\cdot (1 - \\frac{T^m - "
"T_{0}^m}{T_{\\mathrm{melt}}^m - T_{0}^m})$), "
"as shown in Mareau et al. (Mechanics of Materials 143, 2020)"});
}

Expand Down
130 changes: 119 additions & 11 deletions src/mat/4C_mat_inelastic_defgrad_factors.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1643,6 +1643,9 @@ void Mat::InelasticDefgradTransvIsotropElastViscoplast::pre_evaluate(
const Teuchos::ParameterList& params, const EvaluationContext& context, const int gp,
const int eleGID)
{
// save parameter list
params_ = params;

// set Gauss Point
gp_ = gp;

Expand All @@ -1662,7 +1665,7 @@ void Mat::InelasticDefgradTransvIsotropElastViscoplast::pre_evaluate(
// call pre_evaluate method of the time step quantities
time_step_quantities_.pre_evaluate(gp);
// call preevaluate method of the viscoplastic law
viscoplastic_law_->pre_evaluate(gp);
viscoplastic_law_->pre_evaluate(params, gp);
}

/*--------------------------------------------------------------------*
Expand Down Expand Up @@ -2231,7 +2234,7 @@ Mat::InelasticDefgradTransvIsotropElastViscoplast::evaluate_state_quantity_deriv
dNpdMe_sym_dev = Core::LinAlg::Voigt::modify_voigt_representation(temp6x6, 1.0, 2.0);

// compute the relevant derivatives of the plastic strain rate
Core::LinAlg::Matrix<2, 1> evoEqFunctionDers =
InelasticDefgradTransvIsotropElastViscoplastUtils::PlasticStrainRateDerivs evoEqFunctionDers =
viscoplastic_law_->evaluate_derivatives_of_plastic_strain_rate(equiv_stress, plastic_strain,
dt, parameter()->max_plastic_strain_deriv_incr(), err_status);

Expand All @@ -2249,8 +2252,8 @@ Mat::InelasticDefgradTransvIsotropElastViscoplast::evaluate_state_quantity_deriv
}

// compute derivatives of the plastic strain rate
state_quantity_derivatives.curr_dpsr_dequiv_stress = evoEqFunctionDers(0);
state_quantity_derivatives.curr_dpsr_depsp = evoEqFunctionDers(1);
state_quantity_derivatives.curr_dpsr_dequiv_stress = evoEqFunctionDers.deriv_equiv_stress;
state_quantity_derivatives.curr_dpsr_depsp = evoEqFunctionDers.deriv_plastic_strain;


if (eval_type == ViscoplastUtils::StateQuantityDerivEvalType::PlasticStrainRateDerivsOnly)
Expand Down Expand Up @@ -2427,7 +2430,7 @@ void Mat::InelasticDefgradTransvIsotropElastViscoplast::evaluate_additional_cmat
time_step_quantities_.current_plastic_strain[gp_]);

Core::LinAlg::Matrix<10, 10> jacMat(Core::LinAlg::Initialization::zero);
viscoplastic_law_->pre_evaluate(gp_); // set last_substep <- last_
viscoplastic_law_->pre_evaluate(params_, gp_); // set last_substep <- last_
Copy link
Copy Markdown
Contributor

@rjoussen rjoussen Mar 22, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this pre_evaluate needed here?

Copy link
Copy Markdown
Contributor Author

@dragos-ana dragos-ana Mar 27, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good question. The original motivation was to set correct history variables in the case of substepping for the subsequent evaluations. Since my approach has deviated a lot from the original substepping approach, I do not actively use this anymore. But I think that an additional pre_evaluate call does not break anything, although I will have to look at this in more detail at some point. Thanks for spotting it!

jacMat = calculate_jacobian(CredM, current_sol,
time_step_quantities_.last_plastic_defgrad_inverse[gp_],
time_step_quantities_.last_plastic_strain[gp_], time_step_tracker_.dt, err_status);
Expand Down Expand Up @@ -2565,6 +2568,7 @@ void Mat::InelasticDefgradTransvIsotropElastViscoplast::evaluate_inverse_inelast
{
time_step_quantities_.current_plastic_defgrad_inverse[gp_] = iFinM;
time_step_quantities_.current_plastic_strain[gp_] = plastic_strain_pred;
time_step_quantities_.current_equiv_stress[gp_] = state_quantities_.curr_equiv_stress;
time_step_quantities_.current_rightCG[gp_] = CredM;
time_step_quantities_.current_defgrad[gp_] = FredM;
}
Expand Down Expand Up @@ -2596,6 +2600,7 @@ void Mat::InelasticDefgradTransvIsotropElastViscoplast::evaluate_inverse_inelast
{
time_step_quantities_.current_plastic_defgrad_inverse[gp_] = iFinM;
time_step_quantities_.current_plastic_strain[gp_] = sol(9);
time_step_quantities_.current_equiv_stress[gp_] = state_quantities_.curr_equiv_stress;
time_step_quantities_.current_rightCG[gp_] = CredM;
time_step_quantities_.current_defgrad[gp_] = FredM;
}
Expand Down Expand Up @@ -2714,7 +2719,7 @@ Mat::InelasticDefgradTransvIsotropElastViscoplast::calculate_local_newton_loop_r
Core::LinAlg::Matrix<3, 3> resFM(Core::LinAlg::Initialization::zero);
double resepsp = 0.0;

// compute residuals (standard substepping)
// compute residuals (standard time integration)
if (parameter()->timint_type() == ViscoplastUtils::TimIntType::standard)
{
// calculate residual of the equation for inelastic defgrad
Expand All @@ -2733,11 +2738,31 @@ Mat::InelasticDefgradTransvIsotropElastViscoplast::calculate_local_newton_loop_r
last_FinM.invert(last_iFinM);
Core::LinAlg::Matrix<3, 3> T(Core::LinAlg::Initialization::zero);
T.multiply_nn(1.0, last_FinM, iFinM, 0.0);
Core::LinAlg::MatrixFunctErrorType log_err_status{
Core::LinAlg::MatrixFunctErrorType::no_errors};
Core::LinAlg::Matrix<3, 3> logT = Core::LinAlg::matrix_log(T, log_err_status);
FOUR_C_ASSERT_ALWAYS(log_err_status == Core::LinAlg::MatrixFunctErrorType::no_errors,
"Matrix logarithm evaluation failed!");
Core::LinAlg::MatrixFunctErrorType log_err_status =
Core::LinAlg::MatrixFunctErrorType::no_errors;
Core::LinAlg::Matrix<3, 3> logT{Core::LinAlg::Initialization::zero};
if (parameter()->mat_log_calc_method() ==
Core::LinAlg::MatrixLogCalcMethod::inv_scal_square) // evaluation using the inverse
// scaling and squaring
// method
{
// when computing the matrix logarithm with the inverse scaling
// and squaring, we also save the resulting Pade
// order via the dedicated pointer. This will be helpful when we
// compute the derivative - we want consistent Pade orders for the
// evaluations of functions and their derivatives.
logT = Core::LinAlg::matrix_log(
T, log_err_status, matrix_exp_log_utils_.pade_order, parameter()->mat_log_calc_method());
}
else // evaluation using other provided methods
{
logT = Core::LinAlg::matrix_log(T, log_err_status, parameter()->mat_log_calc_method());
}
if (log_err_status != Core::LinAlg::MatrixFunctErrorType::no_errors)
{
err_status = ViscoplastUtils::ErrorType::failed_matrix_log_evaluation;
return Core::LinAlg::Matrix<10, 1>{Core::LinAlg::Initialization::zero};
}

// calculate residual of the equation for inelastic defgrad
resFM.update(1.0, logT, dt, state_quantities_.curr_lpM, 0.0);
Expand Down Expand Up @@ -3313,5 +3338,88 @@ std::string Mat::InelasticDefgradTransvIsotropElastViscoplast::get_error_info(
}


/*--------------------------------------------------------------------*
*--------------------------------------------------------------------*/
void Mat::InelasticDefgradTransvIsotropElastViscoplast::register_output_data_names(
std::unordered_map<std::string, int>& names_and_size) const
{
names_and_size["inverse_plastic_defgrad"] = 9;
names_and_size["plastic_strain"] = 1;
names_and_size["equiv_stress"] = 1;
names_and_size["defgrad"] = 9;
names_and_size["rightCG"] = 9;
viscoplastic_law_->register_output_data_names(names_and_size);
}

/*--------------------------------------------------------------------*
*--------------------------------------------------------------------*/
bool Mat::InelasticDefgradTransvIsotropElastViscoplast::evaluate_output_data(
const std::string& name, Core::LinAlg::SerialDenseMatrix& data) const
{
// auxiliaries
Core::LinAlg::Matrix<9, 1> temp9x1{Core::LinAlg::Initialization::zero};

if (name == "inverse_plastic_defgrad")
{
for (int gp = 0;
gp < static_cast<int>(time_step_quantities_.current_plastic_defgrad_inverse.size()); ++gp)
{
Core::LinAlg::Voigt::matrix_3x3_to_9x1(
time_step_quantities_.current_plastic_defgrad_inverse[gp], temp9x1);

for (int col = 0; col < 9; ++col)
{
data(gp, col) = temp9x1(col);
}
}
return true;
}
else if (name == "plastic_strain")
{
for (int gp = 0; gp < static_cast<int>(time_step_quantities_.current_plastic_strain.size());
++gp)
{
data(gp, 0) = time_step_quantities_.current_plastic_strain[gp];
}
return true;
}
else if (name == "equiv_stress")
{
for (int gp = 0; gp < static_cast<int>(time_step_quantities_.current_equiv_stress.size()); ++gp)
{
data(gp, 0) = time_step_quantities_.current_equiv_stress[gp];
}
return true;
}
else if (name == "defgrad")
{
for (int gp = 0; gp < static_cast<int>(time_step_quantities_.current_defgrad.size()); ++gp)
{
Core::LinAlg::Voigt::matrix_3x3_to_9x1(time_step_quantities_.current_defgrad[gp], temp9x1);

for (int col = 0; col < 9; ++col)
{
data(gp, col) = temp9x1(col);
}
}
return true;
}
else if (name == "rightCG")
{
for (int gp = 0; gp < static_cast<int>(time_step_quantities_.current_rightCG.size()); ++gp)
{
Core::LinAlg::Voigt::matrix_3x3_to_9x1(time_step_quantities_.current_rightCG[gp], temp9x1);


for (int col = 0; col < 9; ++col)
{
data(gp, col) = temp9x1(col);
}
}
return true;
}

return viscoplastic_law_->evaluate_output_data(name, data);
}

FOUR_C_NAMESPACE_CLOSE
32 changes: 32 additions & 0 deletions src/mat/4C_mat_inelastic_defgrad_factors.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -654,6 +654,32 @@ namespace Mat

virtual void unpack_inelastic(Core::Communication::UnpackBuffer& data) = 0;

/*!
* @brief Register names of the internal data that should be saved during runtime output
*
* @param[out] name_and_size Unordered map of names of the data with the respective vector size
*/
virtual void register_output_data_names(
std::unordered_map<std::string, int>& names_and_size) const
{
}

/*!
* @brief Evaluate internal data for every Gauss point saved for output during runtime
* output
*
* @param[in] name Name of the data to export
* @param[out] data NUMGPxNUMDATA Matrix holding the data
*
* @return true if data is set by the material, otherwise false
*/
virtual bool evaluate_output_data(
const std::string& name, Core::LinAlg::SerialDenseMatrix& data) const
{
return false;
}


private:
/// material parameters
Core::Mat::PAR::Parameter* params_;
Expand Down Expand Up @@ -1326,6 +1352,12 @@ namespace Mat

void unpack_inelastic(Core::Communication::UnpackBuffer& buffer) override;

void register_output_data_names(
std::unordered_map<std::string, int>& names_and_size) const override;

bool evaluate_output_data(
const std::string& name, Core::LinAlg::SerialDenseMatrix& data) const override;

/*! @brief Evaluate the current state variables based on a given right Cauchy-Green
* deformation tensor, given inverse plastic deformation gradient and given equivalent
* plastic strain
Expand Down
19 changes: 19 additions & 0 deletions src/mat/4C_mat_inelastic_defgrad_factors_service.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -170,12 +170,17 @@ void Mat::InelasticDefgradTransvIsotropElastViscoplastUtils::TimeStepQuantities:
current_plastic_strain.resize(1, 0.0); // value irrelevant at this point
last_substep_plastic_strain.resize(1, 0.0);

// update last_ and current_ values of the equivalent stress
last_equiv_stress.resize(1, 0.0);
current_equiv_stress.resize(1, 0.0); // value irrelevant at this point

// default values of the right CG tensor: unit tensor
last_rightCG.resize(1, id3x3);
current_rightCG.resize(1, id3x3); // value irrelevant at this point

// default value for the current deformation gradient: zero tensor \f$ \boldsymbol{0} f$ (to make
// sure that the inverse inelastic deformation gradient is evaluated in the first method call)
last_defgrad.resize(1, Core::LinAlg::Matrix<3, 3>{id3x3});
current_defgrad.resize(1, Core::LinAlg::Matrix<3, 3>{Core::LinAlg::Initialization::zero});
}

Expand All @@ -190,6 +195,7 @@ void Mat::InelasticDefgradTransvIsotropElastViscoplastUtils::TimeStepQuantities:
"you attempt to set it to {}",
last_plastic_strain.size(), numgp);


// default values of the inverse plastic deformation gradient for ALL Gauss Points
last_plastic_defgrad_inverse.resize(numgp, last_plastic_defgrad_inverse[0]);
current_plastic_defgrad_inverse.resize(numgp,
Expand All @@ -205,7 +211,12 @@ void Mat::InelasticDefgradTransvIsotropElastViscoplastUtils::TimeStepQuantities:
last_rightCG.resize(numgp, last_rightCG[0]);
current_rightCG.resize(numgp, last_rightCG[0]); // value irrelevant at this point

// default values of the equivalent stress for ALL Gauss Points
last_equiv_stress.resize(numgp, last_equiv_stress[0]);
current_equiv_stress.resize(numgp, current_equiv_stress[0]); // value irrelevant at this point

// default values of the deformation gradient
last_defgrad.resize(numgp, last_defgrad[0]);
current_defgrad.resize(numgp, current_defgrad[0]);
}

Expand All @@ -224,10 +235,12 @@ void Mat::InelasticDefgradTransvIsotropElastViscoplastUtils::TimeStepQuantities:
void Mat::InelasticDefgradTransvIsotropElastViscoplastUtils::TimeStepQuantities::update()
{
// update history variables for the next time step
last_defgrad = current_defgrad;
last_rightCG = current_rightCG;
last_plastic_defgrad_inverse = current_plastic_defgrad_inverse;
last_substep_plastic_defgrad_inverse = current_plastic_defgrad_inverse;
last_plastic_strain = current_plastic_strain;
last_equiv_stress = current_equiv_stress;
last_substep_plastic_strain = current_plastic_strain;
}

Expand All @@ -237,9 +250,11 @@ void Mat::InelasticDefgradTransvIsotropElastViscoplastUtils::TimeStepQuantities:
void Mat::InelasticDefgradTransvIsotropElastViscoplastUtils::TimeStepQuantities::pack(
Core::Communication::PackBuffer& data) const
{
add_to_pack(data, last_defgrad);
add_to_pack(data, last_rightCG);
add_to_pack(data, last_plastic_defgrad_inverse);
add_to_pack(data, last_plastic_strain);
add_to_pack(data, last_equiv_stress);
add_to_pack(data, last_substep_plastic_defgrad_inverse);
add_to_pack(data, last_substep_plastic_strain);
}
Expand All @@ -250,9 +265,11 @@ void Mat::InelasticDefgradTransvIsotropElastViscoplastUtils::TimeStepQuantities:
Core::Communication::UnpackBuffer& buffer)
{
// extract last values
extract_from_pack(buffer, last_defgrad);
extract_from_pack(buffer, last_rightCG);
extract_from_pack(buffer, last_plastic_defgrad_inverse);
extract_from_pack(buffer, last_plastic_strain);
extract_from_pack(buffer, last_equiv_stress);
extract_from_pack(buffer, last_substep_plastic_defgrad_inverse);
extract_from_pack(buffer, last_substep_plastic_strain);

Expand All @@ -263,6 +280,8 @@ void Mat::InelasticDefgradTransvIsotropElastViscoplastUtils::TimeStepQuantities:
last_plastic_defgrad_inverse[0]); // value irrelevant
current_plastic_strain.resize(last_plastic_strain.size(),
last_plastic_strain[0]); // value irrelevant
current_equiv_stress.resize(last_equiv_stress.size(),
last_equiv_stress[0]); // value irrelevant

// set evaluated deformation gradient to 0, to make sure that the inverse inelastic deformation
// gradient is evaluated fully after the restart
Expand Down
Loading
Loading