Hermitian beam elements - discretization and boundary conditions

[davidrudlstorfer]

Question / discussion for all beam experts of 4C (@isteinbrecher, …)

I currently do not know if the following behaviour is a feature or a bug - therefore, I am opening this discussion to summarize and discuss everything

Prerequisites

Based on the improved convergence rate (see https://www.sciencedirect.com/science/article/pii/S0020768317303372) it is preferred to use cubic Hermite polynomials for the centreline discretization and quadratic Lagrange polynomials for the triad interpolation (-> HERM2LINE3 elements in 4C) for some problem types. This type of beam discretization is also used in the following.

The beam discretization looks like

node-1        node-2        node-3
  o-------------o--------------o

with one element being defined with

--------------------------------------STRUCTURE ELEMENTS
1 BEAM3R HERM2LINE3 1 3 2 MAT 1 TRIADS 0.0 0.0 0.0  0.0 0.0 0.0  0.0 0.0 0.0

Problem setup

In my setup I simply want to move a beam element as a rigid body through space. Therefore, I want to fix all DOF's of the beam element. Further, I want to record the resulting forces/moments of the boundary condition with the TAG monitor_reaction

Due to the TAG monitor_reaction only working with --DESIGN POINT DIRICH CONDITIONS and not DESIGN LINE DIRICH CONDITIONS I only use DBCs on individual nodes/points.

Not working / "broken" version

For a long time I am using the following setup which is reproduced in the following minimal working example.

Minimal working example "broken"

--------------------------------------PROBLEM TYPE
PROBLEMTYPE                           Structure
--------------------------------------STRUCTURAL DYNAMIC
LINEAR_SOLVER                         1
INT_STRATEGY                          Standard
DYNAMICTYPE                           OneStepTheta
NLNSOL                                fullnewton
TIMESTEP                              0.1
NUMSTEP                               10
MAXTIME                               1.0
NEGLECTINERTIA                        yes
--------------------------------------STRUCTURAL DYNAMIC/ONESTEPTHETA
THETA                                 1.0
--------------------------------------SOLVER 1
NAME                                  Structure_Solver
SOLVER                                UMFPACK
--------------------------------------IO/MONITOR STRUCTURE DBC
INTERVAL_STEPS                        1
--------------------------------------MATERIALS
MAT 1 MAT_BeamReissnerElastHyper YOUNG 100 POISSONRATIO 0.3 DENS 1.0 CROSSAREA 1.0 SHEARCORR 1.0 MOMINPOL 1.0 MOMIN2 1.0 MOMIN3 1.0
--------------------------------------FUNCT1
SYMBOLIC_FUNCTION_OF_SPACE_TIME a
VARIABLE 0 NAME a TYPE linearinterpolation NUMPOINTS 2 TIMES 0.0 1.0 VALUES 0.0 10.0
--------------------------------------DESIGN POINT DIRICH CONDITIONS
E 1 NUMDOF 9 ONOFF 1 1 1 1 1 1 1 1 1 VAL 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 FUNCT 1 0 0 0 0 0 0 0 0 TAG monitor_reaction
--------------------------------------DNODE-NODE TOPOLOGY
NODE 1 DNODE 1
NODE 2 DNODE 1
NODE 3 DNODE 1
--------------------------------------NODE COORDS
NODE 1 COORD  0.0 0.0 0.0
NODE 2 COORD  5.0 0.0 0.0
NODE 3 COORD 10.0 0.0 0.0
--------------------------------------STRUCTURE ELEMENTS
1 BEAM3R HERM2LINE3 1 3 2 MAT 1 TRIADS 0.0 0.0 0.0  0.0 0.0 0.0  0.0 0.0 0.0

Just recently I found out that this now produces a not wanted / wrong force acting on the beam:

Force output "broken"

step                    time                ref_area               curr_area                     f_x                     f_y                     f_z                     m_x                     m_y                    m_z
 0  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
 1  1.0000000000000001e-01  0.0000000000000000e+00  0.0000000000000000e+00 -2.0512820512823026e+01  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
 2  2.0000000000000001e-01  0.0000000000000000e+00  0.0000000000000000e+00 -4.1025641025646053e+01  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
 3  3.0000000000000004e-01  0.0000000000000000e+00  0.0000000000000000e+00 -6.1538461538469086e+01  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
 4  4.0000000000000002e-01  0.0000000000000000e+00  0.0000000000000000e+00  4.6834570403689554e+01  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
 5  5.0000000000000000e-01  0.0000000000000000e+00  0.0000000000000000e+00  2.6321749890866531e+01  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
 6  5.9999999999999998e-01  0.0000000000000000e+00  0.0000000000000000e+00  5.8089293780435112e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
 7  6.9999999999999996e-01  0.0000000000000000e+00  0.0000000000000000e+00 -1.4703891134779512e+01  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
 8  7.9999999999999993e-01  0.0000000000000000e+00  0.0000000000000000e+00 -3.5216711647602523e+01  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
 9  8.9999999999999991e-01  0.0000000000000000e+00  0.0000000000000000e+00 -5.5729532160425549e+01  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
10  9.9999999999999989e-01  0.0000000000000000e+00  0.0000000000000000e+00  5.2643499781733112e+01  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00

As you can see we now get a random force in the x-direction of the DBC.

Working version

Just recently I found out that this appearing force in x-direction can be attributed to the DBC. With the following version the random force vanishes

Minimal working example "working"

--------------------------------------PROBLEM TYPE
PROBLEMTYPE                           Structure
--------------------------------------STRUCTURAL DYNAMIC
LINEAR_SOLVER                         1
INT_STRATEGY                          Standard
DYNAMICTYPE                           OneStepTheta
NLNSOL                                fullnewton
TIMESTEP                              0.1
NUMSTEP                               10
MAXTIME                               1.0
NEGLECTINERTIA                        yes
--------------------------------------STRUCTURAL DYNAMIC/ONESTEPTHETA
THETA                                 1.0
--------------------------------------SOLVER 1
NAME                                  Structure_Solver
SOLVER                                UMFPACK
--------------------------------------IO/MONITOR STRUCTURE DBC
INTERVAL_STEPS                        1
--------------------------------------MATERIALS
MAT 1 MAT_BeamReissnerElastHyper YOUNG 100 POISSONRATIO 0.3 DENS 1.0 CROSSAREA 1.0 SHEARCORR 1.0 MOMINPOL 1.0 MOMIN2 1.0 MOMIN3 1.0
--------------------------------------FUNCT1
SYMBOLIC_FUNCTION_OF_SPACE_TIME a
VARIABLE 0 NAME a TYPE linearinterpolation NUMPOINTS 2 TIMES 0.0 1.0 VALUES 0.0 10.0
--------------------------------------DESIGN POINT DIRICH CONDITIONS
E 1 NUMDOF 9 ONOFF 1 1 1 1 1 1 1 1 1 VAL 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 FUNCT 1 0 0 0 0 0 0 0 0 TAG monitor_reaction
E 2 NUMDOF 3 ONOFF 1 1 1 VAL 0.0 0.0 0.0 FUNCT 0 0 0
--------------------------------------DNODE-NODE TOPOLOGY
NODE 1 DNODE 1
NODE 2 DNODE 2
NODE 3 DNODE 1
--------------------------------------NODE COORDS
NODE 1 COORD  0.0 0.0 0.0
NODE 2 COORD  5.0 0.0 0.0
NODE 3 COORD 10.0 0.0 0.0
--------------------------------------STRUCTURE ELEMENTS
1 BEAM3R HERM2LINE3 1 3 2 MAT 1 TRIADS 0.0 0.0 0.0  0.0 0.0 0.0  0.0 0.0 0.0

Force output "working"

step                    time                ref_area               curr_area                     f_x                     f_y                     f_z                     m_x                     m_y                    m_z
 0  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
 1  1.0000000000000001e-01  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
 2  2.0000000000000001e-01  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
 3  3.0000000000000004e-01  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
 4  4.0000000000000002e-01  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
 5  5.0000000000000000e-01  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
 6  5.9999999999999998e-01  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
 7  6.9999999999999996e-01  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
 8  7.9999999999999993e-01  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
 9  8.9999999999999991e-01  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00
10  9.9999999999999989e-01  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00

This now produces a force response as expected - every force/moment vanishes.

Differences

The difference between both setup lies within the DBC

The broken version applies the DBC to all three nodes of the beam with:

--------------------------------------DESIGN POINT DIRICH CONDITIONS
E 1 NUMDOF 9 ONOFF 1 1 1 1 1 1 1 1 1 VAL 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 FUNCT 1 0 0 0 0 0 0 0 0 TAG monitor_reaction

The working version applies the DBC seperately to the outer and inner nodes like:

--------------------------------------DESIGN POINT DIRICH CONDITIONS
E 1 NUMDOF 9 ONOFF 1 1 1 1 1 1 1 1 1 VAL 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 FUNCT 1 0 0 0 0 0 0 0 0 TAG monitor_reaction
E 2 NUMDOF 3 ONOFF 1 1 1 VAL 0.0 0.0 0.0 FUNCT 0 0 0

After some further tinkering I found out why the first version produces the seemingly wrong results. If you apply the function to the center node (i.e., rotational degree of freedom) the results are identical to the wrong/"broken" version:

--------------------------------------DESIGN POINT DIRICH CONDITIONS
E 1 NUMDOF 9 ONOFF 1 1 1 1 1 1 1 1 1 VAL 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 FUNCT 1 0 0 0 0 0 0 0 0 TAG monitor_reaction
E 2 NUMDOF 3 ONOFF 1 1 1 VAL 1.0 0.0 0.0 FUNCT 1 0 0 TAG monitor_reaction

Expected behaviour

I always thought that 4C is "smart enough" / that it is implemented so I can apply DBCs to all nodes of a HERM2LINE3 element at once.

I know that the outer nodes have 9 degrees of freedom (displacement, rotation, tangents) and the center nodes have 3 degrees of freedom (rotation).

As you can see above - if you apply the DBC to all nodes at once, 4C applies the DBC of the displacement onto the rotational degrees of freedom of the center node. Therefore, the results are wrong.

Question / Discussion

Is it wanted that one has to split up the DBC for the outer/center nodes of a HERM2LINE3 beam element?

In my opinion it should be possible to use one DBC and 4C correctly applies the boundary conditions to the correct DOFs

[davidrudlstorfer]

And one thing regarding input (@sebproell)

Would it make sense to check if too much DOF's are given in a DBC?

For example I have a node with three degrees of freedom - but I can define a random amount of DOF's for that node

E 2 NUMDOF 7 ONOFF 1 1 1 2 3 4 5 VAL 0.0 0.0 0.0 1.0 2.0 3.0 4.0 FUNCT 0 0 0 1 1 1 1

This would also work for a node with 3 DOFs and the remainder is just ignored

[sebproell]

I think the function that applies the DBCs should for sure check that the number of dofs is exactly what is expected. So, this is only tangentially related to the input mechanism.

[mayrmt]

The problem somehow originates form the fact, that we just specify a list of DBC values. Their assignment to a certain DOF is not part of the DBC input line, but has to be done by the physical field. Ultimately, the interpretation, which entry in the DBC input line acts on which DOF, is not given explicitly, but is only encoded implicitly in each physical field.

For the example with beams described above, the first three entries in the list of DBC values is interpreted as "displacement" values for the end nodes of the element. However, the middle node interprets them as rotational DOFs.

So far, this is just an observation. When thinking about it, it might be better to encode this directly in the DBC input line, however this is a major change and I cannot even "4C" 😉 its implications right now ... we might end up with having to implement physics-specific DBC input lines... not so sure, if this would be a bright idea.

Nevertheless, these are my two cents for this discussion so far.

[isteinbrecher]

That is a very good question. With DBCs in 4C we define which discrete DOF(s) of a node are constrained. The physical meaning of those DOF are encoded in the elements and can not (at least with the standard DBC conditions) directly be constrained.

Take your example of the Beam3rHerm2Line3, if you contain the 1st DOF for all nodes of the element, you will constrain the x-displacement on the corner nodes and the x-component of the spin vector for the middle node. For a KL element the same DBC would mean something completely different. I don’t consider this as a bug, but rather the expected behaviour. I understand that at first glance it might be nice to b able to specify DBC on „physical“ quantities, but that is not as straight forward as one might think. In the beam case, one has to go very deep in the respective beam formulations to get the DBCs right and catch all corner cases (for some theories you simply can’t constrain some quantities). The good news, in practice all of this is no big deal if you follow a few simple rules:
Only apply DBC on the „boundary“ nodes of a beam element. Don’t touch the middle node. If you feel like you have to constrain the middle node, then in 99% of the cases your DBCs on the middle node make no physical sense and you should reconsider your model.

A few comments on your example:

• Be aware that the DBC force monitor outputs forces (and their resulting moment around the origin) for nodes where the first 3 DOFs are displacements in x, y, z. All other node types likely lead to wrong outputs. I assume, the „random“ forces you see are moments enforcing your rotational DBCs in the middle nodes.
• Are you sure that our OST implementation can handle rotational DOFs?
• Why do you expect the force to be 0? You start with a beam at rest and accelerate it to a constant velocity in the first step, so I would expect a non-zero force there.

[davidrudlstorfer]

Thanks for your exhaustive answer @isteinbrecher!

Only apply DBC on the „boundary“ nodes of a beam element. Don’t touch the middle node. If you feel like you have to constrain the middle node, then in 99% of the cases your DBCs on the middle node make no physical sense and you should reconsider your model.

Maybe we can briefly talk about this in person some time soon

Regarding your comments

Be aware that the DBC force monitor outputs forces (and their resulting moment around the origin) for nodes where the first 3 DOFs are displacements in x, y, z. All other node types likely lead to wrong outputs. I assume, the „random“ forces you see are moments enforcing your rotational DBCs in the middle nodes.

Yes that's exactly what was happening

Are you sure that our OST implementation can handle rotational DOFs?

As far as I know - yes. Once I use the necessary problem type Polymer_Network to also activate the Brownian Dynamics Framework only Backwards Euler is supported because the Brownian Dynamics Framework only works with Backwards Euler. You can see this setup in all beam3r_herm2line3_backweuler_browndyn_*** test cases.

Why do you expect the force to be 0? You start with a beam at rest and accelerate it to a constant velocity in the first step, so I would expect a non-zero force there.

I expect the force to be zero because I also activated NEGLECTINERTIA yes. Therefore, no active forces should be on the beam element. The neglect inertia option is necessary because with the Backwards Euler scheme we are not able to determine the accelerations. Note that if some other combination is executed (i.e., other time integration approach or the neglectinertia option is not activated) - the brownian dynamics framework won't work.