Skip to content
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

[StructuralMechanicsApplication] Add initial state and stress calculation to truss element and constitutive law #12525

Merged
merged 11 commits into from
Sep 4, 2024
Merged

[StructuralMechanicsApplication] Add initial state and stress calculation to truss element and constitutive law #12525

Show file tree
Hide file tree
Changes from 10 commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
57d346b
Added the addition of the initial state to the truss constitutive mod…
rfaasse Jul 10, 2024
8788bf0
Reverted some refactoring steps because they give issues with the sco…
rfaasse Jul 10, 2024
56024fe
Added a test for the calculation of the PK2 stress
rfaasse Jul 10, 2024
67dc960
Refactored the additions to the truss element
rfaasse Jul 10, 2024
1ab70da
Added pre-stress, to make it more comparable to the non-linear variant
rfaasse Jul 10, 2024
e4953a3
Cleaning up and making the prestress vector 1 length again
rfaasse Jul 10, 2024
a9d80a0
Changed license for the new test file
rfaasse Jul 10, 2024
a47d09f
Processed review comments and removed some empty statements (redundan…
rfaasse Jul 11, 2024
2a93bc5
Add [[nodiscard]] to a few stub functions in the truss test
rfaasse Jul 11, 2024
7e28fe4
Processing review comments
rfaasse Jul 12, 2024
efd8e60
Merge remote-tracking branch 'origin/master' into geo/12496-add-initi…
rfaasse Sep 3, 2024
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
2 changes: 2 additions & 0 deletions applications/StructuralMechanicsApplication/custom_constitutive/truss_constitutive_law.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -138,9 +138,11 @@ array_1d<double, 3 > & TrussConstitutiveLaw::CalculateValue(
//************************************************************************************
void TrussConstitutiveLaw::CalculateMaterialResponsePK2(Parameters& rValues)
{
AddInitialStrainVectorContribution(rValues.GetStrainVector());
Vector& stress_vector = rValues.GetStressVector();
if (stress_vector.size() != 1) stress_vector.resize(1, false);
stress_vector[0] = this->CalculateStressElastic(rValues);
AddInitialStressVectorContribution(stress_vector);
Copy link
Member

Choose a reason for hiding this comment

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

no initial strains? :)

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Good point, we do not need it for our application at the moment, but added it now for completeness (as well as to the unit test)

}
//************************************************************************************
//************************************************************************************
Expand Down
88 changes: 48 additions & 40 deletions applications/StructuralMechanicsApplication/custom_elements/truss_element_linear_3D2N.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,7 +69,7 @@ TrussElementLinear3D2N::CreateElementStiffnessMatrix(
void TrussElementLinear3D2N::AddPrestressLinear(
VectorType& rRightHandSideVector)
{
KRATOS_TRY;
KRATOS_TRY
BoundedMatrix<double, msLocalSize, msLocalSize> transformation_matrix =
ZeroMatrix(msLocalSize, msLocalSize);
CreateTransformationMatrix(transformation_matrix);
Expand Down Expand Up @@ -145,17 +145,8 @@ void TrussElementLinear3D2N::CalculateOnIntegrationPoints(

array_1d<double, msDimension> temp_internal_stresses = ZeroVector(msDimension);

ConstitutiveLaw::Parameters Values(GetGeometry(),GetProperties(),rCurrentProcessInfo);
Vector temp_strain = ZeroVector(1);
Vector temp_stress = ZeroVector(1);
temp_strain[0] = CalculateLinearStrain();
Values.SetStrainVector(temp_strain);
Values.SetStressVector(temp_stress);
mpConstitutiveLaw->CalculateMaterialResponse(Values,ConstitutiveLaw::StressMeasure_PK2);



truss_forces[0] = (temp_stress[0] + prestress) * A;
const auto stress = CalculateStressFromLinearStrain(rCurrentProcessInfo);
truss_forces[0] = (stress + prestress) * A;

rOutput[0] = truss_forces;
}
Expand All @@ -179,6 +170,15 @@ void TrussElementLinear3D2N::CalculateOnIntegrationPoints(
Strain[2] = 0.00;
rOutput[0] = Strain;
}
else if (rVariable == PK2_STRESS_VECTOR) {
auto stress = CalculateStressFromLinearStrain(rCurrentProcessInfo);
if (GetProperties().Has(TRUSS_PRESTRESS_PK2)) {
stress += GetProperties()[TRUSS_PRESTRESS_PK2];
}

rOutput[0] = ScalarVector(1, stress);
}

KRATOS_CATCH("")
}

Expand All @@ -189,20 +189,22 @@ void TrussElementLinear3D2N::WriteTransformationCoordinates(
BoundedVector<double, TrussElement3D2N::msLocalSize>
& rReferenceCoordinates)
{
KRATOS_TRY;
rReferenceCoordinates = ZeroVector(msLocalSize);
rReferenceCoordinates[0] = GetGeometry()[0].X0();
rReferenceCoordinates[1] = GetGeometry()[0].Y0();
rReferenceCoordinates[2] = GetGeometry()[0].Z0();
rReferenceCoordinates[3] = GetGeometry()[1].X0();
rReferenceCoordinates[4] = GetGeometry()[1].Y0();
rReferenceCoordinates[5] = GetGeometry()[1].Z0();
KRATOS_CATCH("");
KRATOS_TRY
const auto& r_initial_position_node_1 = GetGeometry()[0].GetInitialPosition();
rReferenceCoordinates[0] = r_initial_position_node_1.X();
rReferenceCoordinates[1] = r_initial_position_node_1.Y();
rReferenceCoordinates[2] = r_initial_position_node_1.Z();

const auto& r_initial_position_node_2 = GetGeometry()[1].GetInitialPosition();
rReferenceCoordinates[3] = r_initial_position_node_2.X();
rReferenceCoordinates[4] = r_initial_position_node_2.Y();
rReferenceCoordinates[5] = r_initial_position_node_2.Z();
KRATOS_CATCH("")
}

double TrussElementLinear3D2N::CalculateLinearStrain()
{
KRATOS_TRY;
KRATOS_TRY

Vector current_disp = ZeroVector(msLocalSize);
GetValuesVector(current_disp);
Expand All @@ -215,26 +217,19 @@ double TrussElementLinear3D2N::CalculateLinearStrain()
const double e = (current_disp[3]-current_disp[0])/length_0;

return e;
KRATOS_CATCH("");
KRATOS_CATCH("")
}


void TrussElementLinear3D2N::UpdateInternalForces(BoundedVector<double,msLocalSize>& rInternalForces, const ProcessInfo& rCurrentProcessInfo)
{
KRATOS_TRY;
KRATOS_TRY

Vector temp_internal_stresses = ZeroVector(msLocalSize);
ConstitutiveLaw::Parameters Values(GetGeometry(),GetProperties(),rCurrentProcessInfo);
const auto stress = CalculateStressFromLinearStrain(rCurrentProcessInfo);

Vector temp_strain = ZeroVector(1);
Vector temp_stress = ZeroVector(1);
temp_strain[0] = CalculateLinearStrain();
Values.SetStrainVector(temp_strain);
Values.SetStressVector(temp_stress);
mpConstitutiveLaw->CalculateMaterialResponse(Values,ConstitutiveLaw::StressMeasure_PK2);

temp_internal_stresses[0] = -1.0*temp_stress[0];
temp_internal_stresses[3] = temp_stress[0];
Vector temp_internal_stresses = ZeroVector(msLocalSize);
temp_internal_stresses[0] = -1.0 * stress;
temp_internal_stresses[3] = stress;

rInternalForces = temp_internal_stresses*GetProperties()[CROSS_AREA];

Expand All @@ -245,20 +240,20 @@ void TrussElementLinear3D2N::UpdateInternalForces(BoundedVector<double,msLocalSi

rInternalForces = prod(transformation_matrix, rInternalForces);

KRATOS_CATCH("");
KRATOS_CATCH("")
}

void TrussElementLinear3D2N::FinalizeSolutionStep(const ProcessInfo& rCurrentProcessInfo)
{
KRATOS_TRY;
KRATOS_TRY
ConstitutiveLaw::Parameters Values(GetGeometry(),GetProperties(),rCurrentProcessInfo);
Vector temp_strain = ZeroVector(1);
Vector temp_stress = ZeroVector(1);
temp_strain[0] = CalculateLinearStrain();
Values.SetStrainVector(temp_strain);
Values.SetStressVector(temp_stress);
mpConstitutiveLaw->FinalizeMaterialResponse(Values,ConstitutiveLaw::StressMeasure_PK2);
KRATOS_CATCH("");
KRATOS_CATCH("")
}

double TrussElementLinear3D2N::ReturnTangentModulus1D(const ProcessInfo& rCurrentProcessInfo)
Expand All @@ -268,15 +263,28 @@ double TrussElementLinear3D2N::ReturnTangentModulus1D(const ProcessInfo& rCurren
KRATOS_CATCH("")
}

double TrussElementLinear3D2N::CalculateStressFromLinearStrain(const ProcessInfo &rCurrentProcessInfo)
{
ConstitutiveLaw::Parameters Values(GetGeometry(),GetProperties(),rCurrentProcessInfo);

Vector temp_strain = ScalarVector(1, CalculateLinearStrain());
Vector temp_stress = ZeroVector(1);

Values.SetStrainVector(temp_strain);
Values.SetStressVector(temp_stress);
mpConstitutiveLaw->CalculateMaterialResponse(Values,ConstitutiveLaw::StressMeasure_PK2);

return temp_stress[0];
}

void TrussElementLinear3D2N::save(Serializer& rSerializer) const
{
KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, TrussElement3D2N);
KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, TrussElement3D2N)
rSerializer.save("mConstitutiveLaw", mpConstitutiveLaw);
}
void TrussElementLinear3D2N::load(Serializer& rSerializer)
{
KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, TrussElement3D2N);
KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, TrussElement3D2N)
rSerializer.load("mConstitutiveLaw", mpConstitutiveLaw);
}

Expand Down
1 change: 1 addition & 0 deletions applications/StructuralMechanicsApplication/custom_elements/truss_element_linear_3D2N.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -126,6 +126,7 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) TrussElementLinear3D2N : publ
using TrussElement3D2N::ReturnTangentModulus1D;

private:
double CalculateStressFromLinearStrain(const ProcessInfo &rCurrentProcessInfo);

friend class Serializer;
void save(Serializer& rSerializer) const override;
Expand Down
49 changes: 43 additions & 6 deletions applications/StructuralMechanicsApplication/tests/cpp_tests/test_truss.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,12 +34,12 @@ class StubBilinearLaw : public ConstitutiveLaw

StubBilinearLaw() = default;

ConstitutiveLaw::Pointer Clone() const override
[[nodiscard]] ConstitutiveLaw::Pointer Clone() const override
{
return std::make_shared<StubBilinearLaw>(*this);
}

SizeType GetStrainSize() const override
[[nodiscard]] SizeType GetStrainSize() const override
{
return 1;
}
Expand Down Expand Up @@ -100,9 +100,7 @@ std::shared_ptr<StubBilinearLaw> CreateStubBilinearLaw(double Elongation, double
}


namespace Kratos
{
namespace Testing
namespace Kratos::Testing
{

void AddDisplacementDofsElement(ModelPart& rModelPart){
Expand Down Expand Up @@ -346,5 +344,44 @@ namespace Testing
KRATOS_EXPECT_TRUE(p_constitutive_law->IsInitialized())
}

}
KRATOS_TEST_CASE_IN_SUITE(TrussElementLinear3D2N_CalculatesPK2Stress, KratosStructuralMechanicsFastSuite)
{
Model current_model;
auto &r_model_part = current_model.CreateModelPart("ModelPart",1);
r_model_part.GetProcessInfo().SetValue(DOMAIN_SIZE, 3);
r_model_part.AddNodalSolutionStepVariable(DISPLACEMENT);

// Set the element properties
auto p_elem_prop = r_model_part.CreateNewProperties(0);
constexpr auto youngs_modulus = 2.0e+06;
p_elem_prop->SetValue(YOUNG_MODULUS, youngs_modulus);
const auto &r_clone_cl = KratosComponents<ConstitutiveLaw>::Get("TrussConstitutiveLaw");
p_elem_prop->SetValue(CONSTITUTIVE_LAW, r_clone_cl.Clone());

// Create the test element
constexpr double directional_length = 2.0;
auto p_node_1 = r_model_part.CreateNewNode(1, 0.0, 0.0, 0.0);
auto p_node_2 = r_model_part.CreateNewNode(2, directional_length, directional_length, directional_length);

AddDisplacementDofsElement(r_model_part);

std::vector<ModelPart::IndexType> element_nodes {1,2};
auto p_element = r_model_part.CreateNewElement("TrussLinearElement3D2N", 1, element_nodes, p_elem_prop);
const auto& r_process_info = r_model_part.GetProcessInfo();
p_element->Initialize(r_process_info); // Initialize the element to initialize the constitutive law

constexpr auto induced_strain = 0.1;
p_element->GetGeometry()[1].FastGetSolutionStepValue(DISPLACEMENT) += ScalarVector(3, induced_strain * directional_length);

std::vector<Vector> stress_vector;
p_element->CalculateOnIntegrationPoints(PK2_STRESS_VECTOR, stress_vector, r_process_info);

constexpr auto expected_stress = induced_strain * youngs_modulus;
KRATOS_EXPECT_DOUBLE_EQ(expected_stress, stress_vector[0][0]);

constexpr auto pre_stress = 1.0e5;
p_element->GetProperties().SetValue(TRUSS_PRESTRESS_PK2, pre_stress);
p_element->CalculateOnIntegrationPoints(PK2_STRESS_VECTOR, stress_vector, r_process_info);
KRATOS_EXPECT_DOUBLE_EQ(expected_stress + pre_stress, stress_vector[0][0]);
}
}
65 changes: 65 additions & 0 deletions applications/StructuralMechanicsApplication/tests/cpp_tests/test_truss_constitutive_law.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,65 @@
// KRATOS ___| | | |
// \___ \ __| __| | | __| __| | | __| _` | |
// | | | | | ( | | | | ( | |
// _____/ \__|_| \__,_|\___|\__|\__,_|_| \__,_|_| MECHANICS
//
// License: BSD License
// license: StructuralMechanicsApplication/license.txt
//
// Main authors: Richard Faasse
//

#include "structural_mechanics_fast_suite.h"
#include "custom_constitutive/truss_constitutive_law.h"

namespace Kratos::Testing
{

KRATOS_TEST_CASE_IN_SUITE(TrussConstitutiveLaw_CalculatesLinearElasticStress, KratosStructuralMechanicsFastSuite) {
TrussConstitutiveLaw law;
ConstitutiveLaw::Parameters parameters;
constexpr auto induced_strain = 0.005;
Vector temp_strain = ScalarVector(1, induced_strain);
Vector temp_stress = ZeroVector(1);
parameters.SetStrainVector(temp_strain);
parameters.SetStressVector(temp_stress);

Properties properties;
constexpr auto youngs_modulus = 1.0e6;
properties[YOUNG_MODULUS] = youngs_modulus;
parameters.SetMaterialProperties(properties);
law.CalculateMaterialResponsePK2(parameters);

constexpr double expected_stress = induced_strain * youngs_modulus;
KRATOS_EXPECT_EQ(expected_stress, parameters.GetStressVector()[0]);
}

KRATOS_TEST_CASE_IN_SUITE(TrussConstitutiveLaw_CalculatesLinearElasticStress_WithInitialState, KratosStructuralMechanicsFastSuite) {
TrussConstitutiveLaw law;
constexpr auto induced_strain = 0.005;
Vector temp_strain = ScalarVector(1, induced_strain);
Vector temp_stress = ZeroVector(1);

ConstitutiveLaw::Parameters parameters;
parameters.SetStrainVector(temp_strain);
parameters.SetStressVector(temp_stress);

Properties properties;
constexpr auto youngs_modulus = 1.0e6;
properties[YOUNG_MODULUS] = youngs_modulus;
parameters.SetMaterialProperties(properties);

constexpr auto initial_stress = 5000.0;
constexpr auto initial_strain = 0.01;
const auto p_initial_state = Kratos::make_intrusive<InitialState>(
ScalarVector(1, initial_strain), ScalarVector(1, initial_stress));
law.SetInitialState(p_initial_state);

law.CalculateMaterialResponsePK2(parameters);

constexpr double expected_stress =
(induced_strain - initial_strain) * youngs_modulus + initial_stress;
KRATOS_EXPECT_EQ(expected_stress, parameters.GetStressVector()[0]);
}

}
Loading