From 9a31f244670e232a932c800b84b874f3cb7e54ee Mon Sep 17 00:00:00 2001 From: markelov Date: Tue, 17 Dec 2024 16:05:59 +0100 Subject: [PATCH] - added CalculateOnIntegrationPoints(Vector,.. for LinearTimoshenkoBeamElement2D2N - added TIMOSHENKO_BEAM_PRESTRESS_PK2 - added unit test for this function - moved InitializeConstitutiveLawValuesForStressCalculation to StructuralMechanicsElementUtilities --- .../timoshenko_beam_element_2D2N.cpp | 69 ++++++++------ .../timoshenko_beam_element_2D2N.h | 12 +++ .../truss_elements/linear_truss_element.cpp | 19 +--- .../truss_elements/linear_truss_element.h | 2 - ...tructural_mechanics_python_application.cpp | 1 + ...structural_mechanics_element_utilities.cpp | 18 ++++ .../structural_mechanics_element_utilities.h | 6 ++ .../structural_mechanics_application.cpp | 2 + ...ctural_mechanics_application_variables.cpp | 2 + ...ructural_mechanics_application_variables.h | 1 + .../tests/cpp_tests/test_timoshenko_beams.cpp | 93 +++++++++++++++++++ 11 files changed, 180 insertions(+), 45 deletions(-) create mode 100644 applications/StructuralMechanicsApplication/tests/cpp_tests/test_timoshenko_beams.cpp diff --git a/applications/StructuralMechanicsApplication/custom_elements/beam_elements/timoshenko_beam_element_2D2N.cpp b/applications/StructuralMechanicsApplication/custom_elements/beam_elements/timoshenko_beam_element_2D2N.cpp index 87b221a853de..c6728f151700 100644 --- a/applications/StructuralMechanicsApplication/custom_elements/beam_elements/timoshenko_beam_element_2D2N.cpp +++ b/applications/StructuralMechanicsApplication/custom_elements/beam_elements/timoshenko_beam_element_2D2N.cpp @@ -850,18 +850,11 @@ void LinearTimoshenkoBeamElement2D2N::CalculateOnIntegrationPoints( const auto &r_geometry = GetGeometry(); ConstitutiveLaw::Parameters cl_values(r_geometry, r_props, rProcessInfo); - auto &r_cl_options = cl_values.GetOptions(); - r_cl_options.Set(ConstitutiveLaw::COMPUTE_STRESS , true); - r_cl_options.Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR, false); - const double length = CalculateLength(); const double Phi = StructuralMechanicsElementUtilities::CalculatePhi(r_props, length); - - // Let's initialize the cl values VectorType strain_vector(strain_size), stress_vector(strain_size); - strain_vector.clear(); - cl_values.SetStrainVector(strain_vector); - cl_values.SetStressVector(stress_vector); + StructuralMechanicsElementUtilities::InitializeConstitutiveLawValuesForStressCalculation(cl_values, strain_vector, stress_vector); + VectorType nodal_values(mat_size); GetNodalValuesVector(nodal_values); @@ -880,18 +873,10 @@ void LinearTimoshenkoBeamElement2D2N::CalculateOnIntegrationPoints( const auto &r_geometry = GetGeometry(); ConstitutiveLaw::Parameters cl_values(r_geometry, r_props, rProcessInfo); - auto &r_cl_options = cl_values.GetOptions(); - r_cl_options.Set(ConstitutiveLaw::COMPUTE_STRESS , true); - r_cl_options.Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR, false); - const double length = CalculateLength(); const double Phi = StructuralMechanicsElementUtilities::CalculatePhi(r_props, length); - - // Let's initialize the cl values VectorType strain_vector(strain_size), stress_vector(strain_size); - strain_vector.clear(); - cl_values.SetStrainVector(strain_vector); - cl_values.SetStressVector(stress_vector); + StructuralMechanicsElementUtilities::InitializeConstitutiveLawValuesForStressCalculation(cl_values, strain_vector, stress_vector); VectorType nodal_values(mat_size); GetNodalValuesVector(nodal_values); @@ -910,18 +895,10 @@ void LinearTimoshenkoBeamElement2D2N::CalculateOnIntegrationPoints( const auto &r_geometry = GetGeometry(); ConstitutiveLaw::Parameters cl_values(r_geometry, r_props, rProcessInfo); - auto &r_cl_options = cl_values.GetOptions(); - r_cl_options.Set(ConstitutiveLaw::COMPUTE_STRESS , true); - r_cl_options.Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR, false); - const double length = CalculateLength(); const double Phi = StructuralMechanicsElementUtilities::CalculatePhi(r_props, length); - - // Let's initialize the cl values VectorType strain_vector(strain_size), stress_vector(strain_size); - strain_vector.clear(); - cl_values.SetStrainVector(strain_vector); - cl_values.SetStressVector(stress_vector); + StructuralMechanicsElementUtilities::InitializeConstitutiveLawValuesForStressCalculation(cl_values, strain_vector, stress_vector); VectorType nodal_values(mat_size); GetNodalValuesVector(nodal_values); @@ -964,6 +941,44 @@ void LinearTimoshenkoBeamElement2D2N::CalculateOnIntegrationPoints( } } +void LinearTimoshenkoBeamElement2D2N::CalculateOnIntegrationPoints( + const Variable& rVariable, + std::vector& rOutput, + const ProcessInfo& rProcessInfo + ) +{ + const auto& integration_points = IntegrationPoints(GetIntegrationMethod()); + const SizeType strain_size = mConstitutiveLawVector[0]->GetStrainSize(); + const SizeType mat_size = GetDoFsPerNode() * GetGeometry().size(); + rOutput.resize(integration_points.size()); + const auto &r_props = GetProperties(); + + if (rVariable == PK2_STRESS_VECTOR) { + + const auto &r_geometry = GetGeometry(); + + ConstitutiveLaw::Parameters cl_values(r_geometry, r_props, rProcessInfo); + const double length = CalculateLength(); + const double Phi = StructuralMechanicsElementUtilities::CalculatePhi(r_props, length); + + VectorType strain_vector(strain_size), stress_vector(strain_size); + StructuralMechanicsElementUtilities::InitializeConstitutiveLawValuesForStressCalculation(cl_values, strain_vector, stress_vector); + + VectorType nodal_values(mat_size); + GetNodalValuesVector(nodal_values); + + // Loop over the integration points + for (SizeType integration_point = 0; integration_point < integration_points.size(); ++integration_point) { + CalculateGeneralizedStrainsVector(strain_vector, length, Phi, integration_points[integration_point].X(), nodal_values); + mConstitutiveLawVector[integration_point]->CalculateMaterialResponsePK2(cl_values); + auto stress_vector = cl_values.GetStressVector(); + if ( this->GetProperties().Has(TIMOSHENKO_BEAM_PRESTRESS_PK2)) { + stress_vector += this->GetProperties()[TIMOSHENKO_BEAM_PRESTRESS_PK2]; + } + rOutput[integration_point] = stress_vector; + } + } +} /***********************************************************************************/ /***********************************************************************************/ diff --git a/applications/StructuralMechanicsApplication/custom_elements/beam_elements/timoshenko_beam_element_2D2N.h b/applications/StructuralMechanicsApplication/custom_elements/beam_elements/timoshenko_beam_element_2D2N.h index 599f5d06811e..4f60a257f298 100644 --- a/applications/StructuralMechanicsApplication/custom_elements/beam_elements/timoshenko_beam_element_2D2N.h +++ b/applications/StructuralMechanicsApplication/custom_elements/beam_elements/timoshenko_beam_element_2D2N.h @@ -387,6 +387,18 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) LinearTimoshenkoBeamElement2D const ProcessInfo& rCurrentProcessInfo ) override; + /** + * @brief Calculate a double Variable on the Element Constitutive Law + * @param rVariable The variable we want to get + * @param rOutput The values obtained in the integration points + * @param rCurrentProcessInfo the current process info instance + */ + void CalculateOnIntegrationPoints( + const Variable& rVariable, + std::vector& rOutput, + const ProcessInfo& rCurrentProcessInfo + ) override; + /** * @brief Get on rVariable Constitutive Law from the element * @param rVariable The variable we want to get diff --git a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element.cpp b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element.cpp index 8f8a194cde68..dfab5adf8fe5 100644 --- a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element.cpp +++ b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element.cpp @@ -733,7 +733,7 @@ void LinearTrussElement::CalculateOnIntegrationPoints( ConstitutiveLaw::Parameters cl_values(GetGeometry(), GetProperties(), rProcessInfo); VectorType strain_vector(1), stress_vector(1); MatrixType C(1,1); - InitializeConstitutiveLawValues(cl_values, strain_vector, stress_vector, C); + StructuralMechanicsElementUtilities::InitializeConstitutiveLawValuesForStressCalculation(cl_values, strain_vector, stress_vector, C); const double length = CalculateLength(); SystemSizeBoundedArrayType nodal_values(SystemSize); @@ -756,7 +756,7 @@ void LinearTrussElement::CalculateOnIntegrationPoints( ConstitutiveLaw::Parameters cl_values(GetGeometry(), GetProperties(), rProcessInfo); VectorType strain_vector(1), stress_vector(1); MatrixType C(1,1); - InitializeConstitutiveLawValues(cl_values, strain_vector, stress_vector, C); + StructuralMechanicsElementUtilities::InitializeConstitutiveLawValuesForStressCalculation(cl_values, strain_vector, stress_vector, C); const double length = CalculateLength(); SystemSizeBoundedArrayType nodal_values(SystemSize); @@ -773,19 +773,6 @@ void LinearTrussElement::CalculateOnIntegrationPoints( } } -template -void LinearTrussElement::InitializeConstitutiveLawValues(ConstitutiveLaw::Parameters& rValues, - VectorType& rStrainVector, VectorType& rStressVector, MatrixType rConstitutiveMatrix) -{ - rValues.GetOptions().Set(ConstitutiveLaw::COMPUTE_STRESS , true); - rValues.GetOptions().Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR, false); - - rStrainVector.clear(); - rValues.SetStrainVector(rStrainVector); - rValues.SetStressVector(rStressVector); - rValues.SetConstitutiveMatrix(rConstitutiveMatrix); -} - /***********************************************************************************/ /***********************************************************************************/ @@ -803,7 +790,7 @@ void LinearTrussElement::CalculateOnIntegrationPoints( ConstitutiveLaw::Parameters cl_values(GetGeometry(), GetProperties(), rProcessInfo); VectorType strain_vector(1), stress_vector(1); MatrixType constitutive_matrix(1,1); - InitializeConstitutiveLawValues(cl_values, strain_vector, stress_vector, constitutive_matrix); + StructuralMechanicsElementUtilities::InitializeConstitutiveLawValuesForStressCalculation(cl_values, strain_vector, stress_vector, constitutive_matrix); const double length = CalculateLength(); SystemSizeBoundedArrayType nodal_values(SystemSize); diff --git a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element.h b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element.h index a637edc2617f..5b100a37be6d 100644 --- a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element.h +++ b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element.h @@ -477,8 +477,6 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) LinearTrussElement ///@} ///@name Private Operations ///@{ - void InitializeConstitutiveLawValues(ConstitutiveLaw::Parameters& rValues, - VectorType& rStrainVector, VectorType& rStressVector, MatrixType rConstitutiveMatrix); ///@} ///@name Private Access diff --git a/applications/StructuralMechanicsApplication/custom_python/structural_mechanics_python_application.cpp b/applications/StructuralMechanicsApplication/custom_python/structural_mechanics_python_application.cpp index 2e59f0fddb4f..3fa928e6767a 100644 --- a/applications/StructuralMechanicsApplication/custom_python/structural_mechanics_python_application.cpp +++ b/applications/StructuralMechanicsApplication/custom_python/structural_mechanics_python_application.cpp @@ -91,6 +91,7 @@ PYBIND11_MODULE(KratosStructuralMechanicsApplication,m) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m,I33) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m,IZ) KRATOS_REGISTER_IN_PYTHON_VARIABLE(m,BEAM_INITIAL_STRAIN_VECTOR) + KRATOS_REGISTER_IN_PYTHON_VARIABLE(m,TIMOSHENKO_BEAM_PRESTRESS_PK2) // semi rigid beam variables KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, ROTATIONAL_STIFFNESS_AXIS_2) diff --git a/applications/StructuralMechanicsApplication/custom_utilities/structural_mechanics_element_utilities.cpp b/applications/StructuralMechanicsApplication/custom_utilities/structural_mechanics_element_utilities.cpp index e4e6937c2ca1..e541fc3146d2 100644 --- a/applications/StructuralMechanicsApplication/custom_utilities/structural_mechanics_element_utilities.cpp +++ b/applications/StructuralMechanicsApplication/custom_utilities/structural_mechanics_element_utilities.cpp @@ -635,6 +635,24 @@ double CalculatePhi(const Properties& rProperties, const double L) /***********************************************************************************/ /***********************************************************************************/ +void InitializeConstitutiveLawValuesForStressCalculation(ConstitutiveLaw::Parameters& rValues, + Vector& rStrainVector, Vector& rStressVector, Matrix& rConstitutiveMatrix) +{ + InitializeConstitutiveLawValuesForStressCalculation(rValues, rStrainVector, rStressVector); + rValues.SetConstitutiveMatrix(rConstitutiveMatrix); +} + +void InitializeConstitutiveLawValuesForStressCalculation(ConstitutiveLaw::Parameters& rValues, + Vector& rStrainVector, Vector& rStressVector) +{ + rValues.GetOptions().Set(ConstitutiveLaw::COMPUTE_STRESS , true); + rValues.GetOptions().Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR, false); + + rStrainVector.clear(); + rValues.SetStrainVector(rStrainVector); + rValues.SetStressVector(rStressVector); +} + } // namespace StructuralMechanicsElementUtilities. } // namespace Kratos. diff --git a/applications/StructuralMechanicsApplication/custom_utilities/structural_mechanics_element_utilities.h b/applications/StructuralMechanicsApplication/custom_utilities/structural_mechanics_element_utilities.h index 0705076958bf..ac68dafd11ed 100644 --- a/applications/StructuralMechanicsApplication/custom_utilities/structural_mechanics_element_utilities.h +++ b/applications/StructuralMechanicsApplication/custom_utilities/structural_mechanics_element_utilities.h @@ -350,5 +350,11 @@ double GetReferenceRotationAngle2D3NBeam(const GeometryType &rGeometry); */ KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) double CalculatePhi(const Properties& rProperties, const double L); +void InitializeConstitutiveLawValuesForStressCalculation(ConstitutiveLaw::Parameters& rValues, + Vector& rStrainVector, Vector& rStressVector, Matrix& rConstitutiveMatrix); + +void InitializeConstitutiveLawValuesForStressCalculation(ConstitutiveLaw::Parameters& rValues, + Vector& rStrainVector, Vector& rStressVector); + } // namespace StructuralMechanicsElementUtilities. } // namespace Kratos. diff --git a/applications/StructuralMechanicsApplication/structural_mechanics_application.cpp b/applications/StructuralMechanicsApplication/structural_mechanics_application.cpp index ea5e0292d16e..107540fcb982 100644 --- a/applications/StructuralMechanicsApplication/structural_mechanics_application.cpp +++ b/applications/StructuralMechanicsApplication/structural_mechanics_application.cpp @@ -347,6 +347,8 @@ void KratosStructuralMechanicsApplication::Register() { KRATOS_REGISTER_VARIABLE(I33) KRATOS_REGISTER_VARIABLE(LUMPED_MASS_ROTATION_COEFFICIENT) KRATOS_REGISTER_VARIABLE(BEAM_INITIAL_STRAIN_VECTOR) + KRATOS_REGISTER_VARIABLE(TIMOSHENKO_BEAM_PRESTRESS_PK2) + // semi rigid beam variables KRATOS_REGISTER_VARIABLE(ROTATIONAL_STIFFNESS_AXIS_2) diff --git a/applications/StructuralMechanicsApplication/structural_mechanics_application_variables.cpp b/applications/StructuralMechanicsApplication/structural_mechanics_application_variables.cpp index 262279e7e901..fc67e0d8b434 100644 --- a/applications/StructuralMechanicsApplication/structural_mechanics_application_variables.cpp +++ b/applications/StructuralMechanicsApplication/structural_mechanics_application_variables.cpp @@ -71,6 +71,8 @@ KRATOS_CREATE_VARIABLE(double, I22) KRATOS_CREATE_VARIABLE(double, I33) KRATOS_CREATE_VARIABLE(double, LUMPED_MASS_ROTATION_COEFFICIENT) KRATOS_CREATE_VARIABLE(Vector, BEAM_INITIAL_STRAIN_VECTOR) +KRATOS_CREATE_VARIABLE(Vector, TIMOSHENKO_BEAM_PRESTRESS_PK2) + // semi rigid beam variables KRATOS_CREATE_VARIABLE(double, ROTATIONAL_STIFFNESS_AXIS_2) diff --git a/applications/StructuralMechanicsApplication/structural_mechanics_application_variables.h b/applications/StructuralMechanicsApplication/structural_mechanics_application_variables.h index 93c0d2f739dd..ed0ab5cfb803 100644 --- a/applications/StructuralMechanicsApplication/structural_mechanics_application_variables.h +++ b/applications/StructuralMechanicsApplication/structural_mechanics_application_variables.h @@ -92,6 +92,7 @@ namespace Kratos KRATOS_DEFINE_APPLICATION_VARIABLE( STRUCTURAL_MECHANICS_APPLICATION,double, I33) KRATOS_DEFINE_APPLICATION_VARIABLE( STRUCTURAL_MECHANICS_APPLICATION,double, LUMPED_MASS_ROTATION_COEFFICIENT) KRATOS_DEFINE_APPLICATION_VARIABLE( STRUCTURAL_MECHANICS_APPLICATION,Vector, BEAM_INITIAL_STRAIN_VECTOR) + KRATOS_DEFINE_APPLICATION_VARIABLE( STRUCTURAL_MECHANICS_APPLICATION,Vector, TIMOSHENKO_BEAM_PRESTRESS_PK2) // Semi rigid beam variables diff --git a/applications/StructuralMechanicsApplication/tests/cpp_tests/test_timoshenko_beams.cpp b/applications/StructuralMechanicsApplication/tests/cpp_tests/test_timoshenko_beams.cpp new file mode 100644 index 000000000000..83dc5e317124 --- /dev/null +++ b/applications/StructuralMechanicsApplication/tests/cpp_tests/test_timoshenko_beams.cpp @@ -0,0 +1,93 @@ +// KRATOS ___| | | | +// \___ \ __| __| | | __| __| | | __| _` | | +// | | | | | ( | | | | ( | | +// _____/ \__|_| \__,_|\___|\__|\__,_|_| \__,_|_| MECHANICS +// +// License: BSD License +// license: StructuralMechanicsApplication/license.txt +// +// Main authors: Gennady Markelov +// + +// Project includes +#include "containers/model.h" +#include "structural_mechanics_fast_suite.h" +#include "structural_mechanics_application_variables.h" + +#include +#include + +namespace Kratos::Testing +{ + +void CreateBeamModel2N(std::string TimoshenkoBeamElementName) + { + 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); + r_model_part.AddNodalSolutionStepVariable(ROTATION_Z); + + // 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); + constexpr auto cross_area = 1.0; + p_elem_prop->SetValue(CROSS_AREA, cross_area); + constexpr auto i33 = 1.0; + p_elem_prop->SetValue(I33, i33); + constexpr auto area_effective_y = 1.0; + p_elem_prop->SetValue(AREA_EFFECTIVE_Y, area_effective_y); + + const auto &r_clone_cl = KratosComponents::Get("TimoshenkoBeamElasticConstitutiveLaw"); + 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); + + for (auto& r_node : r_model_part.Nodes()){ + r_node.AddDof(DISPLACEMENT_X); + r_node.AddDof(DISPLACEMENT_Y); + r_node.AddDof(DISPLACEMENT_Z); + r_node.AddDof(ROTATION_Z); + } + + std::vector element_nodes {1,2}; + auto p_element = r_model_part.CreateNewElement(std::move(TimoshenkoBeamElementName), 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); + p_element->GetGeometry()[1].FastGetSolutionStepValue(ROTATION_Z) += 0.1; + + std::vector stress_vectors; + p_element->CalculateOnIntegrationPoints(PK2_STRESS_VECTOR, stress_vectors, r_process_info); + + constexpr auto expected_stress = induced_strain * youngs_modulus; + constexpr auto expected_shear_stress = -37500.0; + + constexpr auto tolerance = 1.0e-5; + Vector expected_stress_vector(3); + expected_stress_vector <<= expected_stress, 29631.5,expected_shear_stress; + KRATOS_EXPECT_VECTOR_RELATIVE_NEAR(expected_stress_vector, stress_vectors[0], tolerance); + expected_stress_vector <<= expected_stress, 70710.7,expected_shear_stress; + KRATOS_EXPECT_VECTOR_RELATIVE_NEAR(expected_stress_vector, stress_vectors[1], tolerance); + expected_stress_vector <<= expected_stress, 111790,expected_shear_stress; + KRATOS_EXPECT_VECTOR_RELATIVE_NEAR(expected_stress_vector, stress_vectors[2], tolerance); + + Vector pre_stress(3); + pre_stress <<= 1.0e5, 1.0e4, 1.0e3; + p_element->GetProperties().SetValue(TIMOSHENKO_BEAM_PRESTRESS_PK2, pre_stress); + p_element->CalculateOnIntegrationPoints(PK2_STRESS_VECTOR, stress_vectors, r_process_info); + expected_stress_vector += pre_stress; + KRATOS_EXPECT_VECTOR_RELATIVE_NEAR(expected_stress_vector, stress_vectors[2], tolerance); + } + + KRATOS_TEST_CASE_IN_SUITE(LinearTimoshenkoBeam2D2N_CalculatesPK2Stress, KratosStructuralMechanicsFastSuite) + { + CreateBeamModel2N("LinearTimoshenkoBeamElement2D2N"); + } +}