Topic: Nonlinear beams and condensing dofs.

I've tried to make beam2d and beam3d possibly support nonlinear behavior. The way condensing dofs works makes this very hard.

I'd like to replace dof condensing in beam2d and beam3d (the only place where it is used) with slave nodes. As far as I understand it, the goal is to obtain a moment free joint, which is easily achieved with a slave node that locks the displacements, placed on top of an existing node. I think this is much easier to understand what is going on, but mainly, it will enable nonlinear support.



---------------


Also, Beam2D and Beam3D use a hardcoded stiffness matrix, and doens't use the B-matrix ever.
I've tried to replace the computations for the internal forces by integrating in the form "B^T * D * B"
The manual claims it to be a Timoshenko beam, and the header file claims "cubic lateral displacement, quadratic rotations"

For example Beam2D. With B defined as

    B.at(1, 1) =  -1. / l;
    B.at(1, 4) =   1. / l;
    B.at(2, 2) =   ( 6. - 12. * ksi ) / ( l * l * ( 1. + 2. * kappa ) );
    B.at(2, 3) =   ( -2. * ( 2. + kappa ) + 6. * ksi ) / ( l * ( 1. + 2. * kappa ) );
    B.at(2, 5) =   ( -6. + 12. * ksi ) / ( l * l * ( 1. + 2. * kappa ) );
    B.at(2, 6) =   ( -2. * ( 1. - kappa ) + 6. * ksi ) / ( l * ( 1. + 2. * kappa ) );
    B.at(3, 2) =   ( -2. * kappa ) / ( l * ( 1. + 2. * kappa ) );
    B.at(3, 3) =   kappa / ( l * ( 1. + 2. * kappa ) );
    B.at(3, 5) =   2. * kappa / ( l * ( 1. + 2. * kappa ) );
    B.at(3, 6) =   kappa / ( l * ( 1. + 2. * kappa ) );

then   integral  B^T * D * B dV    will not produce the same stiffness matrix as the hardcoded one.
Perhaps some trick to avoid shearlocking? Reverse-engineering this is very hard.
If kappa is close to zero, then this difference is barely noticeable.


For Beam3D i get some bigger trouble though. Again, if kappa is close to zero, then I mostly get the same matrix, except for a few components (the same value in a few places, changes sign):

Tangent = FloatMatrix with dimensions : 12 12
 3.938e+06   0.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00  -3.938e+06   0.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00  
 0.000e+00   4.823e+05   0.000e+00   0.000e+00   0.000e+00   2.412e+05   0.000e+00  -4.823e+05   0.000e+00   0.000e+00   0.000e+00   2.412e+05  
 0.000e+00   0.000e+00   7.973e+05   0.000e+00  -3.987e+05   0.000e+00   0.000e+00   0.000e+00  -7.973e+05   0.000e+00  -3.987e+05   0.000e+00  
 0.000e+00   0.000e+00   0.000e+00   3.920e+04   0.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00  -3.920e+04   0.000e+00   0.000e+00  
 0.000e+00   0.000e+00  -3.987e+05   0.000e+00   2.658e+05   0.000e+00   0.000e+00   0.000e+00   3.987e+05   0.000e+00   1.329e+05   0.000e+00  
 0.000e+00   2.412e+05   0.000e+00   0.000e+00   0.000e+00   1.608e+05   0.000e+00  -2.412e+05   0.000e+00   0.000e+00   0.000e+00   8.039e+04  
-3.938e+06   0.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00   3.938e+06   0.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00  
 0.000e+00  -4.823e+05   0.000e+00   0.000e+00   0.000e+00  -2.412e+05   0.000e+00   4.823e+05   0.000e+00   0.000e+00   0.000e+00  -2.412e+05  
 0.000e+00   0.000e+00  -7.973e+05   0.000e+00   3.987e+05   0.000e+00   0.000e+00   0.000e+00   7.973e+05   0.000e+00   3.987e+05   0.000e+00  
 0.000e+00   0.000e+00   0.000e+00  -3.920e+04   0.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00   3.920e+04   0.000e+00   0.000e+00  
 0.000e+00   0.000e+00  -3.987e+05   0.000e+00   1.329e+05   0.000e+00   0.000e+00   0.000e+00   3.987e+05   0.000e+00   2.658e+05   0.000e+00  
 0.000e+00   2.412e+05   0.000e+00   0.000e+00   0.000e+00   8.039e+04   0.000e+00  -2.412e+05   0.000e+00   0.000e+00   0.000e+00   1.608e+05  
B^T * D * B = FloatMatrix with dimensions : 12 12
 3.937e+06   0.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00  -3.937e+06   0.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00  
 0.000e+00   4.823e+05   0.000e+00   0.000e+00   0.000e+00  -2.412e+05   0.000e+00  -4.823e+05   0.000e+00   0.000e+00   0.000e+00  -2.412e+05  <----- Sign changed!
 0.000e+00   0.000e+00   7.973e+05   0.000e+00  -3.987e+05   0.000e+00   0.000e+00   0.000e+00  -7.973e+05   0.000e+00  -3.987e+05   0.000e+00  
 0.000e+00   0.000e+00   0.000e+00   3.920e+04   0.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00  -3.920e+04   0.000e+00   0.000e+00  
 0.000e+00   0.000e+00  -3.987e+05   0.000e+00   2.658e+05   0.000e+00   0.000e+00   0.000e+00   3.987e+05   0.000e+00   1.329e+05   0.000e+00  
 0.000e+00  -2.412e+05   0.000e+00   0.000e+00   0.000e+00   1.608e+05   0.000e+00   2.412e+05   0.000e+00   0.000e+00   0.000e+00   8.039e+04  
-3.937e+06   0.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00   3.937e+06   0.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00  
 0.000e+00  -4.823e+05   0.000e+00   0.000e+00   0.000e+00   2.412e+05   0.000e+00   4.823e+05   0.000e+00   0.000e+00   0.000e+00   2.412e+05  
 0.000e+00   0.000e+00  -7.973e+05   0.000e+00   3.987e+05   0.000e+00   0.000e+00   0.000e+00   7.973e+05   0.000e+00   3.987e+05   0.000e+00  
 0.000e+00   0.000e+00   0.000e+00  -3.920e+04   0.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00   3.920e+04   0.000e+00   0.000e+00  
 0.000e+00   0.000e+00  -3.987e+05   0.000e+00   1.329e+05   0.000e+00   0.000e+00   0.000e+00   3.987e+05   0.000e+00   2.658e+05   0.000e+00  
 0.000e+00  -2.412e+05   0.000e+00   0.000e+00   0.000e+00   8.039e+04   0.000e+00   2.412e+05   0.000e+00   0.000e+00   0.000e+00   1.608e+05

Which must mean a serious bug in either the B matrix or in the stiffness matrix.

Re: Nonlinear beams and condensing dofs.

Well,  I discovered 2 sign errors in the N matrix, and 2 different sign errors in the B-matrix.
Fixing those, now the B matrix is actually the derivative of N, and K = integral B^T * D * B dx
It all looks correct to me now, and the tests all still pass, so I'm pushing those changes.

3

Re: Nonlinear beams and condensing dofs.

MIkael,

I completely understand your arguments and the problems the static condensation causes in case of  nonlinear functions.
However, I still would prefer to have this traditional linear beam: - simple for (at least our) students, as they use them in some courses and can see that they implement exactly what they know (this is especially good for MSc students, they don't know anything about FEM, they know just displacement method), the condensation makes the model preparation simpler.
I think that with clear documentation, we can live with linear and nonlinear beam elements, in fact there are already many different beam elements: libeam family (mindlin beams with linear interpolation, and their nonlinear variants (Simo -Qu Voc elements).

Re: Nonlinear beams and condensing dofs.

I wanted to handle the thermal expansion in the cross-section (so that all structural elements would support it), which is the initial reason why I made any modifications to the beam elements (since they did it manually themselves). Layered cross-sections also now supports thermal expansion due to these changes.

---------------


Making a "momentless-joint"-node would be as simple as something like this:

...
node 2 coords 3 2.4 0.  0. 
node 3 coords 3 3.8 0.  0. 
slavenode 4 coords 3 3.8 0. 0.  masterdofman 1 3 doftype 3 2 2 0   # Momentless joint
node 5 coords 3 5.8 0.  1.5
...
Beam2d 2 nodes 2 2 3 # Connected to one side of the 3-4 joint
Beam2d 3 nodes 2 4 5 # Connected to the other side of the 3-4 joint

compared to the dof-condensing method:

...
node 2 coords 3 2.4 0.  0. 
node 3 coords 3 3.8 0.  0. 
node 4 coords 3 5.8 0.  1.5
...
Beam2d 2 nodes 2 2 3 DofsToCondense 1 6
Beam2d 3 nodes 2 3 4 DofsToCondense 1 3 

Surely this is easier for everyone?
It seems much more logical to me, to model a joint by connecting two nodes.
One could even implement a "jointnode" that just connects only displacements which would just turn it into the trivial

...
node 2 coords 3 2.4 0.  0. 
node 3 coords 3 3.8 0.  0. 
jointnode 4 coords 3 3.8 0. 0.  master 3
node 5 coords 3 5.8 0.  1.5
...
Beam2d 2 nodes 2 2 3 # Connected to one side of the 3-4 joint
Beam2d 3 nodes 2 4 5 # Connected to the other side of the 3-4 joint

But most importantly, in the output, you actually do get the right values when you use the joint-node. With condensation, you can get the rotation of only one of the elements connected to the node.
What happens with the beam forces and displacements I can't even tell. Perhaps they are correct? I don't know. After messing about with the internal forces and tangent stiffness it's hard to tell what is happening.

Most importantly, this method will give you the actual correct moments in the beam, and the correct rotations/displacements for each beam.


-----

I'm perfectly fine with structural elements being hardcoded for linear response (as I mentioned earlier, i was mostly trying to fix bugs in the thermal expansion in layered cross-sections), but I just can't understand the rationale behind dof condensation when a joint is so much easier.

But even for linear element, we would probably want to very the thickness/moment of inertia, so integrating to obtain the stiffness matrix is still much more preferably.

Post's attachments

beam2d_1.in 2.3 kb, 5 downloads since 2014-01-29 

You don't have the permssions to download the attachments of this post.