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.

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] semi-rigid hinge at 3d2n linear beam #12290

Open
wants to merge 13 commits into
base: master
Choose a base branch
from

Conversation

aronnoordam
Copy link
Member

@aronnoordam aronnoordam commented Apr 18, 2024

added semi rigid hinge to linear 3d2n beam element

according to:

[simões 1996] https://www.sciencedirect.com/science/article/abs/pii/0045794995004270,

@philbucher
Copy link
Member

Can you please briefly explain what this does? I dont have access to the paper

@aronnoordam
Copy link
Member Author

@philbucher using the matrix as described in the paper, a hinge with a rotational stiffness can be added to a beam element. where 0 stiffness corresponds to a normal hinge; infinite stiffness corresponds to a normal connection (as it was before).

using the rotational stiffnesses fixity factors for both sides of the beam are calculated ranging from 0 to 1. The resulting matrix (for euler beams) will be:

image

where alpha1 and alpha2 are the fixity factors

@aronnoordam
Copy link
Member Author

fyi, I need this in order to simulate rail joints which have e.g. 90% of the stiffness of the rail itself

@aronnoordam aronnoordam marked this pull request as ready for review April 24, 2024 11:20
@aronnoordam aronnoordam requested a review from a team as a code owner April 24, 2024 11:20
"I33" : 1e-5,
"SEMI_RIGID_ROTATIONAL_STIFFNESS_VECTOR_AXIS_2": [0e0],
"SEMI_RIGID_ROTATIONAL_STIFFNESS_VECTOR_AXIS_3": [1e7],
"SEMI_RIGID_NODE_IDS": [6]
Copy link
Member

Choose a reason for hiding this comment

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

I am pretty hesitant using IDs as input because they are reordered in MPI

How about using Nodal-Data to specify this?

Copy link
Member Author

Choose a reason for hiding this comment

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

that makes sense. I've changed it

{
KRATOS_TRY

if (this->GetProperties().Has(SEMI_RIGID_NODE_IDS)) {
Copy link
Member

Choose a reason for hiding this comment

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

in a restart, should this be done every time?

I guess not no?

Copy link
Member Author

Choose a reason for hiding this comment

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

this cannot be done anymore, as the check for semi-rigidity cannot be done in Initialize anymore, since there is the possibility that the rotational stiffness is not assigned yet (depending on the order of the proces executions

Comment on lines 60 to 67
KRATOS_ERROR_IF_NOT(this->GetProperties().Has(SEMI_RIGID_ROTATIONAL_STIFFNESS_VECTOR_AXIS_2)) <<
"SEMI_RIGID_ROTATIONAL_STIFFNESS_VECTOR_AXIS_2 must be set when SEMI_RIGID_NODE_IDS is used\n";
KRATOS_ERROR_IF_NOT(this->GetProperties().Has(SEMI_RIGID_ROTATIONAL_STIFFNESS_VECTOR_AXIS_3)) <<
"SEMI_RIGID_ROTATIONAL_STIFFNESS_VECTOR_AXIS_3 must be set when SEMI_RIGID_NODE_IDS is used\n";
KRATOS_ERROR_IF_NOT(this->GetProperties()[SEMI_RIGID_NODE_IDS].size() == this->GetProperties()[SEMI_RIGID_ROTATIONAL_STIFFNESS_VECTOR_AXIS_2].size()) <<
"SEMI_RIGID_ROTATIONAL_STIFFNESS_VECTOR_AXIS_2 must have equal size as SEMI_RIGID_NODE_IDS\n";
KRATOS_ERROR_IF_NOT(this->GetProperties()[SEMI_RIGID_NODE_IDS].size() == this->GetProperties()[SEMI_RIGID_ROTATIONAL_STIFFNESS_VECTOR_AXIS_3].size()) <<
"SEMI_RIGID_ROTATIONAL_STIFFNESS_VECTOR_AXIS_3 must have equal size as SEMI_RIGID_NODE_IDS\n";
Copy link
Member

Choose a reason for hiding this comment

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

I think it would make sense to move this to the Check function

"SEMI_RIGID_ROTATIONAL_STIFFNESS_VECTOR_AXIS_3 must have equal size as SEMI_RIGID_NODE_IDS\n";

// create integer list
Vector semi_rigid_node_id_input = this->GetProperties()[SEMI_RIGID_NODE_IDS];
Copy link
Member

Choose a reason for hiding this comment

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

Here is an example of using the nodal data for specifying the input:

const array_1d<double, 3>& nodal_stiffness = this->GetValue(NODAL_DISPLACEMENT_STIFFNESS);

This way you can be sure that the correct nodes have the input
Additionally, for large number of hinges it will it overcrowd your materials.json

Copy link
Member

@philbucher philbucher left a comment

Choose a reason for hiding this comment

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

minor remaining comments, otherwise good to me 👍

"Parameters" : {
"mesh_id" : 0,
"model_part_name" : "Structure.Parts_semi_rigid_hinge",
"variable_name" : "ROTATIONAL_STIFFNESS_AXIS_3",
Copy link
Member

Choose a reason for hiding this comment

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

FYI you can directly write this in the mdpa, up to you

Begin NodalData BOUNDARY
1
2
973
974
End NodalData

Copy link
Member Author

Choose a reason for hiding this comment

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

ah thanks

Comment on lines +72 to +73
KRATOS_CREATE_VARIABLE(double, ROTATIONAL_STIFFNESS_AXIS_2)
KRATOS_CREATE_VARIABLE(double, ROTATIONAL_STIFFNESS_AXIS_3)
Copy link
Member

Choose a reason for hiding this comment

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

just to confirm, those are actual stiffnesses and not bools to enable/disable the hinge?
And they are nodal quantities, not elemental?

Copy link
Member Author

Choose a reason for hiding this comment

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

yes thats correct, its the rotational nodal stiffness in the local direction of the beam

@@ -107,6 +107,9 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) CrBeamElementLinear3D2N : pub

private:

void CalculateRigidityReductionMatrix(
Copy link
Member

Choose a reason for hiding this comment

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

Please add a brief description what this method does

// 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;
Copy link
Member

Choose a reason for hiding this comment

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

if those must be larger than 0 if defined, then this should be checked in Check. Then you can simplify the logic here

Copy link
Member Author

Choose a reason for hiding this comment

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

it is also allowed to be 0, this means that a normal hinge is simulated. In that case, this if statement is still required, else there is a division by 0. I will add a check to check if its lower than 0, because this is not allowed

Comment on lines +338 to +341
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;
}
Copy link
Member

Choose a reason for hiding this comment

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

you can skip the Has check. If you make sure that they have a valid input, and the Geometry from which you get the nodes from is const, then you can simplify to
const double alpha1_zz = 1 / (1 + 3 * E * Iy / (r_geometry[0].GetValue(ROTATIONAL_STIFFNESS_AXIS_2) * L))

Copy link
Member Author

Choose a reason for hiding this comment

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

In order to not have a change of behaviour of the existing beam, the default beam should be without hinge. This means that the default ROTATIONAL_STIFFNESS_AXIS_2 should be infinite. In my oppinion it is cleaner to do a check if the ROTATIONAL_STIFFNESS_AXIS_2 is actually defined, thus keeping the Has.

philbucher
philbucher previously approved these changes Apr 29, 2024
Copy link
Member

@philbucher philbucher left a comment

Choose a reason for hiding this comment

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

Nice 👍

@aronnoordam
Copy link
Member Author

Nice 👍

Thanks :D, however I did not commit the Check and the brief discription yet, I've done that now in the latest commit

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

2 participants