[StructuralMechanicsApplication] semi-rigid hinge at 3d2n linear beam

applications/StructuralMechanicsApplication/custom_elements/cr_beam_element_linear_3D2N.cpp

@@ -51,6 +51,30 @@ CrBeamElementLinear3D2N::Create(IndexType NewId,
 
 CrBeamElementLinear3D2N::~CrBeamElementLinear3D2N() {}
 
+
+int CrBeamElementLinear3D2N::Check(const ProcessInfo& rCurrentProcessInfo) const {
+
+    KRATOS_TRY
+
+    const GeometryType& r_geom = GetGeometry();
+    // verify that the rotational stiffness around axis 2 and 3 are greater or equal than 0.0 if they are given.
+    for (unsigned int i = 0; i < r_geom.size(); ++i) {
+
+        if (r_geom[i].Has(ROTATIONAL_STIFFNESS_AXIS_2)) {
+            KRATOS_ERROR_IF(r_geom[i].GetValue(ROTATIONAL_STIFFNESS_AXIS_2) < 0.0) 
+                << " ROTATIONAL_STIFFNESS_AXIS_2 should be greater or equal than 0.0" << std::endl;
+        }
+        if (r_geom[i].Has(ROTATIONAL_STIFFNESS_AXIS_3)) {
+            KRATOS_ERROR_IF(r_geom[i].GetValue(ROTATIONAL_STIFFNESS_AXIS_3) < 0.0) 
+                << " ROTATIONAL_STIFFNESS_AXIS_3 should be greater or equal than 0.0" << std::endl;
+        }
+    }
+
+    return CrBeamElement3D2N::Check(rCurrentProcessInfo);
+
+    KRATOS_CATCH("")
+}
+
 void CrBeamElementLinear3D2N::CalculateLocalSystem(
     MatrixType& rLeftHandSideMatrix, VectorType& rRightHandSideVector,
     const ProcessInfo& rCurrentProcessInfo)
@@ -98,6 +122,19 @@ void CrBeamElementLinear3D2N::CalculateLeftHandSide(
     rLeftHandSideMatrix = ZeroMatrix(msElementSize, msElementSize);
     noalias(rLeftHandSideMatrix) += CreateElementStiffnessMatrix_Material();
 
+
+    const GeometryType& r_geometry = GetGeometry();
+
+    // reduce rigidity if required
+    for (IndexType i = 0; i < r_geometry.size(); ++i) {
+        if (r_geometry[i].Has(ROTATIONAL_STIFFNESS_AXIS_2) || r_geometry[i].Has(ROTATIONAL_STIFFNESS_AXIS_3)) {
+            BoundedMatrix<double, msElementSize, msElementSize> rigidity_reduction_matrix = ZeroMatrix(msElementSize, msElementSize);
+            CalculateRigidityReductionMatrix(rigidity_reduction_matrix);
+            rLeftHandSideMatrix = element_prod(rLeftHandSideMatrix, rigidity_reduction_matrix);
+            break;
+        }
+    }
+
     //// start static condensation
     if (Has(CONDENSED_DOF_LIST)) {
         Vector dof_list_input = GetValue(CONDENSED_DOF_LIST);
@@ -297,6 +334,114 @@ void CrBeamElementLinear3D2N::CalculateOnIntegrationPoints(
     KRATOS_CATCH("");
 }
 
+
+void CrBeamElementLinear3D2N::CalculateRigidityReductionMatrix(
+    BoundedMatrix<double, msElementSize, msElementSize>& rRigidityReductionMatrix) const
+{
+    KRATOS_TRY
+
+
+    rRigidityReductionMatrix = ZeroMatrix(msElementSize, msElementSize);
+    const double E = GetProperties()[YOUNG_MODULUS];
+    const double L = StructuralMechanicsElementUtilities::CalculateReferenceLength3D2N(*this);
+
+
+    const double Iy = GetProperties()[I22];
+    const double Iz = GetProperties()[I33];
+
+    const GeometryType& r_geometry = GetGeometry();
+
+    // inititalise fixity factors as completely fixed
+    double alpha1_yy = 1;
+    double alpha2_yy = 1;
+
+    double alpha1_zz = 1;
+    double alpha2_zz = 1;
+
+    // set fixity factors, 0 is no fixity, i.e. a hinge; 1 is completely fixed, the resulting stiffness matrix will remain unchanged
+    if (r_geometry[0].Has(ROTATIONAL_STIFFNESS_AXIS_2)) {
+        const double rotational_stiffness_y_1 = r_geometry[0].GetValue(ROTATIONAL_STIFFNESS_AXIS_2);
+        alpha1_zz = (rotational_stiffness_y_1 > 0) ? 1 / (1 + 3 * E * Iy / (rotational_stiffness_y_1 * L)) : 0.0;
+    }
+    if (r_geometry[1].Has(ROTATIONAL_STIFFNESS_AXIS_2)) {
+        const double rotational_stiffness_y_2 = r_geometry[1].GetValue(ROTATIONAL_STIFFNESS_AXIS_2);
+        alpha2_zz = (rotational_stiffness_y_2 > 0) ? 1 / (1 + 3 * E * Iy / (rotational_stiffness_y_2 * L)) : 0.0;
+    }
+    if (r_geometry[0].Has(ROTATIONAL_STIFFNESS_AXIS_3)) {
+        const double rotational_stiffness_z_1 = r_geometry[0].GetValue(ROTATIONAL_STIFFNESS_AXIS_3);
+        alpha1_yy = (rotational_stiffness_z_1 > 0) ? 1 / (1 + 3 * E * Iz / (rotational_stiffness_z_1 * L)) : 0.0;
+
+    }
+    if (r_geometry[1].Has(ROTATIONAL_STIFFNESS_AXIS_3)) {
+        const double rotational_stiffness_z_2 = r_geometry[1].GetValue(ROTATIONAL_STIFFNESS_AXIS_3);
+        alpha2_yy = (rotational_stiffness_z_2 > 0) ? 1 / (1 + 3 * E * Iz / (rotational_stiffness_z_2 * L)) : 0.0;
+    }
+
+    const double yy_denominator = (4 - alpha1_yy * alpha2_yy);
+    const double zz_denominator = (4 - alpha1_zz * alpha2_zz);
+
+    // normal
+    rRigidityReductionMatrix(0, 0) = 1;
+    rRigidityReductionMatrix(0, 6) = 1;
+    rRigidityReductionMatrix(6, 6) = 1;
+    rRigidityReductionMatrix(6, 0) = 1;
+
+    // yy
+    rRigidityReductionMatrix(1, 1) = (alpha1_yy + alpha2_yy + alpha1_yy * alpha2_yy) / yy_denominator;
+    rRigidityReductionMatrix(7, 7) = rRigidityReductionMatrix(1, 1);
+    rRigidityReductionMatrix(1, 7) = rRigidityReductionMatrix(1, 1);
+    rRigidityReductionMatrix(7, 1) = rRigidityReductionMatrix(1, 1);
+
+    // zz
+    rRigidityReductionMatrix(2, 2) = (alpha1_zz + alpha2_zz + alpha1_zz * alpha2_zz) / zz_denominator;
+    rRigidityReductionMatrix(8, 8) = rRigidityReductionMatrix(2, 2);
+    rRigidityReductionMatrix(2, 8) = rRigidityReductionMatrix(2, 2);
+    rRigidityReductionMatrix(8, 2) = rRigidityReductionMatrix(2, 2);
+
+    // tortion
+    rRigidityReductionMatrix(3, 3) = 1;
+    rRigidityReductionMatrix(9, 9) = 1;
+
+    // yy 1,2 - rot z 1
+    rRigidityReductionMatrix(1, 5) = (2 * alpha1_yy + alpha1_yy * alpha2_yy) / yy_denominator;
+    rRigidityReductionMatrix(5, 1) = rRigidityReductionMatrix(1, 5);
+    rRigidityReductionMatrix(5, 7) = rRigidityReductionMatrix(1, 5);
+    rRigidityReductionMatrix(7, 5) = rRigidityReductionMatrix(1, 5);
+
+    // zz 1,2 - rot y 1
+    rRigidityReductionMatrix(2, 4) = (2 * alpha1_zz + alpha1_zz * alpha2_zz) / zz_denominator;
+    rRigidityReductionMatrix(4, 2) = rRigidityReductionMatrix(2, 4);
+    rRigidityReductionMatrix(4, 8) = rRigidityReductionMatrix(2, 4);
+    rRigidityReductionMatrix(8, 4) = rRigidityReductionMatrix(2, 4);
+
+    // yy 1,2 - rot z 2
+    rRigidityReductionMatrix(1, 11) = (2 * alpha2_yy + alpha1_yy * alpha2_yy) / yy_denominator;
+    rRigidityReductionMatrix(11, 1) = rRigidityReductionMatrix(1, 11);
+    rRigidityReductionMatrix(11, 7) = rRigidityReductionMatrix(1, 11);
+    rRigidityReductionMatrix(7, 11) = rRigidityReductionMatrix(1, 11);
+
+    // zz 1,2 - rot y 2
+    rRigidityReductionMatrix(2, 10) = (2 * alpha2_zz + alpha1_zz * alpha2_zz) / zz_denominator;
+    rRigidityReductionMatrix(10, 2) = rRigidityReductionMatrix(2, 10);
+    rRigidityReductionMatrix(10, 8) = rRigidityReductionMatrix(2, 10);
+    rRigidityReductionMatrix(8, 10) = rRigidityReductionMatrix(2, 10);
+
+    // rot z 1,2 
+    rRigidityReductionMatrix(5, 5) = 3 * alpha1_yy / yy_denominator;
+    rRigidityReductionMatrix(11, 11) = 3 * alpha2_yy / yy_denominator;
+    rRigidityReductionMatrix(5, 11) = 3 * alpha1_yy * alpha2_yy / yy_denominator;
+    rRigidityReductionMatrix(11, 5) = rRigidityReductionMatrix(5, 11);
+
+    // rot y 1, 2
+    rRigidityReductionMatrix(4, 4) = 3 * alpha1_zz / zz_denominator;
+    rRigidityReductionMatrix(10, 10) = 3 * alpha2_zz / zz_denominator;
+    rRigidityReductionMatrix(4, 10) = 3 * alpha1_zz * alpha2_zz / zz_denominator;
+    rRigidityReductionMatrix(10, 4) = rRigidityReductionMatrix(4, 10);
+
+    KRATOS_CATCH("");
+
+}
+
 void CrBeamElementLinear3D2N::save(Serializer& rSerializer) const
 {
     KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, CrBeamElement3D2N);

applications/StructuralMechanicsApplication/custom_elements/cr_beam_element_linear_3D2N.hpp

@@ -70,6 +70,8 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) CrBeamElementLinear3D2N : pub
         PropertiesType::Pointer pProperties
     ) const override;
 
+    int Check(const ProcessInfo& rCurrentProcessInfo) const override;
+
     void CalculateLocalSystem(
         MatrixType& rLeftHandSideMatrix,
         VectorType& rRightHandSideVector,
@@ -107,6 +109,12 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) CrBeamElementLinear3D2N : pub
 
 private:
 
+    /**
+    * @brief This function calculates a matrix which is used to reduce rigidity at the beam connections
+    */
+    void CalculateRigidityReductionMatrix(
+        BoundedMatrix<double, msElementSize, msElementSize>& rReductionMatrix) const;
+
     friend class Serializer;
     void save(Serializer& rSerializer) const override;
     void load(Serializer& rSerializer) override;

applications/StructuralMechanicsApplication/custom_python/structural_mechanics_python_application.cpp

@@ -92,6 +92,10 @@ PYBIND11_MODULE(KratosStructuralMechanicsApplication,m)
     KRATOS_REGISTER_IN_PYTHON_VARIABLE(m,IZ)
     KRATOS_REGISTER_IN_PYTHON_VARIABLE(m,BEAM_INITIAL_STRAIN_VECTOR)
 
+    // semi rigid beam variables
+    KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, ROTATIONAL_STIFFNESS_AXIS_2)
+    KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, ROTATIONAL_STIFFNESS_AXIS_3)
+
     // shell generalized variables
     KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, STENBERG_SHEAR_STABILIZATION_SUITABLE )
     KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, SHELL_OFFSET )

applications/StructuralMechanicsApplication/structural_mechanics_application.cpp

@@ -325,6 +325,10 @@ void KratosStructuralMechanicsApplication::Register() {
     KRATOS_REGISTER_VARIABLE(LUMPED_MASS_ROTATION_COEFFICIENT)
     KRATOS_REGISTER_VARIABLE(BEAM_INITIAL_STRAIN_VECTOR)
 
+    // semi rigid beam variables
+    KRATOS_REGISTER_VARIABLE(ROTATIONAL_STIFFNESS_AXIS_2)
+    KRATOS_REGISTER_VARIABLE(ROTATIONAL_STIFFNESS_AXIS_3)
+
     //  Shell generalized variables
     KRATOS_REGISTER_VARIABLE(STENBERG_SHEAR_STABILIZATION_SUITABLE)
     KRATOS_REGISTER_VARIABLE(SHELL_OFFSET)

applications/StructuralMechanicsApplication/structural_mechanics_application_variables.cpp

@@ -68,6 +68,10 @@ KRATOS_CREATE_VARIABLE(double, I33)
 KRATOS_CREATE_VARIABLE(double, LUMPED_MASS_ROTATION_COEFFICIENT)
 KRATOS_CREATE_VARIABLE(Vector, BEAM_INITIAL_STRAIN_VECTOR)
 
+// semi rigid beam variables
+KRATOS_CREATE_VARIABLE(double, ROTATIONAL_STIFFNESS_AXIS_2)
+KRATOS_CREATE_VARIABLE(double, ROTATIONAL_STIFFNESS_AXIS_3)
+
 // Shell generalized variables
 KRATOS_CREATE_VARIABLE(bool, STENBERG_SHEAR_STABILIZATION_SUITABLE)
 KRATOS_CREATE_VARIABLE(double, SHELL_OFFSET)

applications/StructuralMechanicsApplication/structural_mechanics_application_variables.h

@@ -89,6 +89,11 @@ namespace Kratos
     KRATOS_DEFINE_APPLICATION_VARIABLE( STRUCTURAL_MECHANICS_APPLICATION,double, LUMPED_MASS_ROTATION_COEFFICIENT)
     KRATOS_DEFINE_APPLICATION_VARIABLE( STRUCTURAL_MECHANICS_APPLICATION,Vector, BEAM_INITIAL_STRAIN_VECTOR)
 
+
+    // Semi rigid beam variables
+    KRATOS_DEFINE_APPLICATION_VARIABLE( STRUCTURAL_MECHANICS_APPLICATION, double, ROTATIONAL_STIFFNESS_AXIS_2)
+    KRATOS_DEFINE_APPLICATION_VARIABLE( STRUCTURAL_MECHANICS_APPLICATION, double, ROTATIONAL_STIFFNESS_AXIS_3)
+
     // Shell generalized variables
     KRATOS_DEFINE_APPLICATION_VARIABLE( STRUCTURAL_MECHANICS_APPLICATION, bool, STENBERG_SHEAR_STABILIZATION_SUITABLE )
     KRATOS_DEFINE_APPLICATION_VARIABLE( STRUCTURAL_MECHANICS_APPLICATION, double, SHELL_OFFSET)

applications/StructuralMechanicsApplication/tests/beam_test/semi_rigid_linear_3D2NBeamCr_test.mdpa

@@ -0,0 +1,140 @@
+
+Begin Properties 0
+End Properties
+
+Begin Nodes
+   1        0.00000        0.00000        0.00000
+   2        0.1000        0.00000        0.00000
+   3        0.200        0.00000        0.00000
+   4        0.30000        0.00000        0.00000
+   5        0.4000        0.00000        0.00000
+   6        0.500        0.00000        0.00000
+   7        0.60000        0.00000        0.00000
+   8        0.7000        0.00000        0.00000
+   9        0.800        0.00000        0.00000   
+   10       0.900        0.00000        0.00000   
+   11       1.000        0.00000        0.00000   
+End Nodes
+
+
+Begin Elements CrLinearBeamElement3D2N// GUI group identifier: structure
+         1          0         1 			2
+         2          0         2 			3
+         3          0         3 			4
+         4          0         4 			5
+         5          0         5 			6
+         6          0         6 			7
+         7          0         7 			8
+         8          0         8 			9
+         9          0         9 			10
+         10         0        10 			11
+End Elements
+
+Begin Conditions PointLoadCondition3D1N// GUI group identifier: neumann
+1 0 6
+End Conditions
+
+
+Begin SubModelPart Parts_semi_rigid_hinge
+    Begin SubModelPartNodes
+       6
+    End SubModelPartNodes
+    Begin SubModelPartElements
+    End SubModelPartElements
+    Begin SubModelPartConditions
+    End SubModelPartConditions
+End SubModelPart
+Begin SubModelPart Parts_structure // Group structure // Subtree Parts
+    Begin SubModelPartNodes
+        1
+       2
+       3
+       4
+       5
+       6
+       7
+       8
+       9
+       10
+       11
+    End SubModelPartNodes
+    Begin SubModelPartElements
+         1
+         2
+         3
+         4
+         5
+         6
+         7
+         8
+         9
+         10
+    End SubModelPartElements
+    Begin SubModelPartConditions
+    End SubModelPartConditions
+End SubModelPart
+Begin SubModelPart DISPLACEMENT_dirichletXYZ // Group dirichletXYZ // Subtree DISPLACEMENT
+    Begin SubModelPartNodes
+        1
+		11
+    End SubModelPartNodes
+    Begin SubModelPartElements
+    End SubModelPartElements
+    Begin SubModelPartConditions
+    End SubModelPartConditions
+End SubModelPart
+Begin SubModelPart ROTATION_dirrot // Group dirrot // Subtree ROTATION
+    Begin SubModelPartNodes
+    1
+	11
+    End SubModelPartNodes
+    Begin SubModelPartElements
+    End SubModelPartElements
+    Begin SubModelPartConditions
+    End SubModelPartConditions
+End SubModelPart
+Begin SubModelPart ROTATION_dirrotX // Group dirrot // Subtree ROTATION
+    Begin SubModelPartNodes
+    End SubModelPartNodes
+    Begin SubModelPartElements
+    End SubModelPartElements
+    Begin SubModelPartConditions
+    End SubModelPartConditions
+End SubModelPart
+Begin SubModelPart DISPLACEMENT_dirichletYZ // Group dirichletYZ // Subtree DISPLACEMENT
+    Begin SubModelPartNodes
+    End SubModelPartNodes
+    Begin SubModelPartElements
+    End SubModelPartElements
+    Begin SubModelPartConditions
+    End SubModelPartConditions
+End SubModelPart
+Begin SubModelPart DISPLACEMENT_allnodes // Group allnodes // Subtree DISPLACEMENT
+    Begin SubModelPartNodes
+        1
+       2
+       3
+       4
+       5
+       6
+       7
+       8
+       9
+       10
+       11
+    End SubModelPartNodes
+    Begin SubModelPartElements
+    End SubModelPartElements
+    Begin SubModelPartConditions
+    End SubModelPartConditions
+End SubModelPart
+Begin SubModelPart PointLoad3D_neumann // Group neumann // Subtree PointLoad3D
+    Begin SubModelPartNodes
+       6
+    End SubModelPartNodes
+    Begin SubModelPartElements
+    End SubModelPartElements
+    Begin SubModelPartConditions
+         1
+    End SubModelPartConditions
+End SubModelPart

applications/StructuralMechanicsApplication/tests/beam_test/semi_rigid_linear_3D2NBeamCr_test_material.json

@@ -0,0 +1,22 @@
+{
+    "properties" : [{
+        "model_part_name" : "Structure.Parts_structure",
+        "properties_id"   : 1,
+        "Material"        : {
+            "constitutive_law" : {
+                "name" : "BeamConstitutiveLaw"
+            },
+            "Variables"        : {
+                "DENSITY"           : 7850.0,
+                "YOUNG_MODULUS"     : 210000000000.0,
+                "POISSON_RATIO"     : 0.30,
+                "CROSS_AREA"        : 0.00001,
+                "TORSIONAL_INERTIA" : 1e-5,
+                "I22"               : 1e-5,
+                "I33"               : 1e-5
+            },
+            "Tables"           : {}
+        }
+    }]
+}
+