Skip to content

Commit

Permalink
Merge pull request #11570 from KratosMultiphysics/guo-viscous-calcula…
Browse files Browse the repository at this point in the history 
…tion-prestressed
  • Loading branch information
@AlejandroCornejo
AlejandroCornejo authored Sep 18, 2023
2 parents b10d5cb + 85c2c88 commit 5a16407
Show file tree
Hide file tree
Showing 3 changed files with 18 additions and 74 deletions.
13 changes: 6 additions & 7 deletions ...eLawsApplication/custom_constitutive/small_strains/viscous/viscous_generalized_kelvin.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -121,8 +121,8 @@ void ViscousGeneralizedKelvin<TElasticBehaviourLaw>::ComputeViscoElasticity(Cons
const Flags& r_flags = rValues.GetOptions();

// We compute the strain in case not provided
Vector& r_strain_vector = rValues.GetStrainVector();
if (r_flags.IsNot( ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN)) {
Vector& r_strain_vector = rValues.GetStrainVector();
this->CalculateValue(rValues, STRAIN, r_strain_vector);
}

Expand All @@ -146,7 +146,6 @@ void ViscousGeneralizedKelvin<TElasticBehaviourLaw>::ComputeViscoElasticity(Cons
Vector elastic_strain(VoigtSize);

// Apply increments
const Vector& r_strain_vector = rValues.GetStrainVector();
for (IndexType i = 0; i < number_sub_increments; ++i) {
noalias(aux) = (std::exp(-dt / delay_time) * prod(inverse_constitutive_matrix, aux_stress_vector)) / delay_time;
noalias(inelastic_strain) = std::exp(-dt / delay_time) * inelastic_strain + aux;
Expand All @@ -158,7 +157,8 @@ void ViscousGeneralizedKelvin<TElasticBehaviourLaw>::ComputeViscoElasticity(Cons
noalias(integrated_stress_vector) = aux_stress_vector;

if (r_flags.Is(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR)) {
rValues.SetConstitutiveMatrix(constitutive_matrix);
Matrix& r_constitutive_matrix = rValues.GetConstitutiveMatrix();
noalias(r_constitutive_matrix) = constitutive_matrix;
}
} else {
// Elastic Matrix
Expand Down Expand Up @@ -210,8 +210,8 @@ void ViscousGeneralizedKelvin<TElasticBehaviourLaw>::FinalizeMaterialResponseCau
const Flags& r_flags = rValues.GetOptions();

// We compute the strain in case not provided
Vector& r_strain_vector = rValues.GetStrainVector();
if (r_flags.IsNot( ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN)) {
Vector& r_strain_vector = rValues.GetStrainVector();
this->CalculateValue(rValues, STRAIN, r_strain_vector);
}

Expand All @@ -233,16 +233,15 @@ void ViscousGeneralizedKelvin<TElasticBehaviourLaw>::FinalizeMaterialResponseCau
Vector elastic_strain(VoigtSize);

// Apply increments
const Vector& r_strain_vector = rValues.GetStrainVector();
for (IndexType i = 0; i < number_sub_increments; ++i) {
noalias(aux) = (std::exp(-dt / delay_time) * prod(inverse_constitutive_matrix, aux_stress_vector)) / delay_time;
noalias(inelastic_strain) = std::exp(-dt / delay_time) * inelastic_strain + aux;
noalias(elastic_strain) = r_strain_vector - inelastic_strain;
noalias(aux_stress_vector) = prod(constitutive_matrix, elastic_strain);
}

mPrevInelasticStrainVector = inelastic_strain;
mPrevStressVector = aux_stress_vector;
noalias(mPrevInelasticStrainVector) = inelastic_strain;
noalias(mPrevStressVector) = aux_stress_vector;
}

/***********************************************************************************/
Expand Down
52 changes: 12 additions & 40 deletions ...LawsApplication/custom_constitutive/small_strains/viscous/viscous_generalized_maxwell.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,32 +101,18 @@ void ViscousGeneralizedMaxwell<TElasticBehaviourLaw>::CalculateMaterialResponseC
/***********************************************************************************/
/***********************************************************************************/

template <class TElasticBehaviourLaw>
void ViscousGeneralizedMaxwell<TElasticBehaviourLaw>::FinalizeSolutionStep(
const Properties &rMaterialProperties,
const GeometryType &rElementGeometry,
const Vector& rShapeFunctionsValues,
const ProcessInfo &rCurrentProcessInfo)
{
// Deprecated
}

/***********************************************************************************/
/***********************************************************************************/

template <class TElasticBehaviourLaw>
void ViscousGeneralizedMaxwell<TElasticBehaviourLaw>::ComputeViscoElasticity(ConstitutiveLaw::Parameters& rValues)
{
const Properties& material_props = rValues.GetMaterialProperties();
Vector& integrated_stress_vector = rValues.GetStressVector(); // To be updated
const ProcessInfo& process_info = rValues.GetProcessInfo();
const double time_step = process_info[DELTA_TIME];
const double time_step = rValues.GetProcessInfo()[DELTA_TIME];
const Flags& r_flags = rValues.GetOptions();

// We compute the strain in case not provided
Vector& r_strain_vector = rValues.GetStrainVector();
if (r_flags.IsNot( ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN)) {
Vector& r_strain_vector = rValues.GetStrainVector();
this->CalculateValue(rValues, STRAIN, r_strain_vector);
BaseType::CalculateValue(rValues, STRAIN, r_strain_vector);
}

// We compute the stress
Expand All @@ -136,9 +122,8 @@ void ViscousGeneralizedMaxwell<TElasticBehaviourLaw>::ComputeViscoElasticity(Con

// Elastic Matrix
Matrix constitutive_matrix;
this->CalculateValue(rValues, CONSTITUTIVE_MATRIX, constitutive_matrix);
BaseType::CalculateValue(rValues, CONSTITUTIVE_MATRIX, constitutive_matrix);

const Vector& r_strain_vector = rValues.GetStrainVector();
const Vector& r_previous_strain = this->GetPreviousStrainVector();
const Vector& r_previous_stress = this->GetPreviousStressVector();
const Vector& r_strain_increment = r_strain_vector - r_previous_strain;
Expand All @@ -149,13 +134,14 @@ void ViscousGeneralizedMaxwell<TElasticBehaviourLaw>::ComputeViscoElasticity(Con
noalias(integrated_stress_vector) = r_previous_stress * std::exp(-time_step / delay_time) + prod(constitutive_matrix, r_auxiliary_strain);

if (r_flags.Is(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR)) {
rValues.SetConstitutiveMatrix(constitutive_matrix);
Matrix& r_constitutive_matrix = rValues.GetConstitutiveMatrix();
noalias(r_constitutive_matrix) = constitutive_matrix;
}
} else {
// Elastic Matrix
if( r_flags.Is( ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR ) ) {
Matrix& r_constitutive_matrix = rValues.GetConstitutiveMatrix();
this->CalculateValue(rValues, CONSTITUTIVE_MATRIX, r_constitutive_matrix);
BaseType::CalculateValue(rValues, CONSTITUTIVE_MATRIX, r_constitutive_matrix);
}
}
}
Expand Down Expand Up @@ -203,9 +189,9 @@ void ViscousGeneralizedMaxwell<TElasticBehaviourLaw>::FinalizeMaterialResponseCa
const Flags& r_flags = rValues.GetOptions();

// We compute the strain in case not provided
Vector& r_strain_vector = rValues.GetStrainVector();
if (r_flags.IsNot( ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN)) {
Vector& r_strain_vector = rValues.GetStrainVector();
this->CalculateValue(rValues, STRAIN, r_strain_vector);
BaseType::CalculateValue(rValues, STRAIN, r_strain_vector);
}


Expand All @@ -214,9 +200,8 @@ void ViscousGeneralizedMaxwell<TElasticBehaviourLaw>::FinalizeMaterialResponseCa

// Elastic Matrix
Matrix constitutive_matrix;
this->CalculateValue(rValues, CONSTITUTIVE_MATRIX, constitutive_matrix);
BaseType::CalculateValue(rValues, CONSTITUTIVE_MATRIX, constitutive_matrix);

const Vector& r_strain_vector = rValues.GetStrainVector();
const Vector& r_previous_strain = this->GetPreviousStrainVector();
const Vector& r_previous_stress = this->GetPreviousStressVector();
const Vector& r_strain_increment = r_strain_vector - r_previous_strain;
Expand All @@ -226,21 +211,8 @@ void ViscousGeneralizedMaxwell<TElasticBehaviourLaw>::FinalizeMaterialResponseCa

noalias(integrated_stress_vector) = r_previous_stress * std::exp(-time_step / delay_time) + prod(constitutive_matrix, r_auxiliary_strain);

mPrevStressVector = integrated_stress_vector;
mPrevStrainVector = r_strain_vector;
}

/***********************************************************************************/
/***********************************************************************************/

template <class TElasticBehaviourLaw>
Vector& ViscousGeneralizedMaxwell<TElasticBehaviourLaw>::CalculateValue(
ConstitutiveLaw::Parameters& rParameterValues,
const Variable<Vector>& rThisVariable,
Vector& rValue
)
{
return BaseType::CalculateValue(rParameterValues, rThisVariable, rValue);
noalias(mPrevStressVector) = integrated_stress_vector;
noalias(mPrevStrainVector) = r_strain_vector;
}

/***********************************************************************************/
Expand Down
27 changes: 0 additions & 27 deletions ...veLawsApplication/custom_constitutive/small_strains/viscous/viscous_generalized_maxwell.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -150,20 +150,6 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) ViscousGeneralizedMaxwell
*/
void CalculateMaterialResponseCauchy(ConstitutiveLaw::Parameters& rValues) override;

/**
* to be called at the end of each solution step
* (e.g. from Element::FinalizeSolutionStep)
* @param rMaterialProperties the Properties instance of the current element
* @param rElementGeometry the geometry of the current element
* @param rShapeFunctionsValues the shape functions values in the current integration point
* @param the current ProcessInfo instance
*/
void FinalizeSolutionStep(
const Properties &rMaterialProperties,
const GeometryType &rElementGeometry,
const Vector &rShapeFunctionsValues,
const ProcessInfo &rCurrentProcessInfo) override;

/**
* Finalize the material response in terms of 1st Piola-Kirchhoff stresses
* @see Parameters
Expand All @@ -188,19 +174,6 @@ class KRATOS_API(CONSTITUTIVE_LAWS_APPLICATION) ViscousGeneralizedMaxwell
*/
void FinalizeMaterialResponseCauchy(ConstitutiveLaw::Parameters& rValues) override;

/**
* @brief Calculates the value of a specified variable (Vector)
* @param rParameterValues the needed parameters for the CL calculation
* @param rThisVariable the variable to be returned
* @param rValue a reference to the returned value
* @param rValue output: the value of the specified variable
*/
Vector& CalculateValue(
ConstitutiveLaw::Parameters& rParameterValues,
const Variable<Vector>& rThisVariable,
Vector& rValue
) override;

/**
* @brief Calculates the value of a specified variable (Matrix)
* @param rParameterValues the needed parameters for the CL calculation
Expand Down

0 comments on commit 5a16407

Please sign in to comment.