Dear all,
I noticed something strange in the results of Beam3d elements. In the attached example I analyzed a cantilever beam with square cross section and 2 loads on the free end in directions -globalZ (-20) and globalY (10):

1 - file two-st.in uses refAngle = 0
2 - file two-st-refNode uses node 3 (coords 0 0 1 - same as globalZ) as reference node for the cross section; it gives the same results as in case 1
3 - file two-st-refNode uses node 3 (coords 0 0 1 - same as globalY) as reference node for the cross section; results looks strange in the DofManager output (I would expect to see the same displacements...)

The issues I found are:
a. Myy (5th and 11th values) has not the correct sign, I would expect the opposite of the one I get;
b. the ouput files are not reporting in any line the correct element number! I read beam element 1 and in the input deck it is defined as number 5
c. is the square cross section definition correct? I read this page but I suspect I'm missing something in it.

Finally, I would definitely prefer to use refAngle instead of refNode, but I'm still confused with the results I get. I'm running OOFEM updated to the last code from Mikael's Git.
Please don't mind the unuseful lines like the 1st cross section or the doubled nodalload.

thanks,
Giovanni

Thanks Mikael.

I was thinking mee too the lines:

    answer.at(5, 3) =   ( 6. - 12. * ksi ) / ( l * l * c1y );
...
    answer.at(6, 2) =  -( 6. - 12. * ksi ) / ( l * l * c1z ); // new

and they seems to have the opposite sign, as can be seen  here.

In the meanwhile, I tested the same files with beamShearCoefficient 1.e30 (as I found in the example files) and with the correct factor (5/6 for the rectangular section).

In this last case, the equilibrium is not respected, and also the local displacements are not correct. I attach some hand calculations to clarify; this is very annoying and I cannot understand where the bug is located (maybe in the same 2 lines above).

Thanks, I'll give it a try.

Yes, I posted the stiffness matrix for 2D beam, we need the B matrix in the 3d case, sorry.

So, the problem has nothing to do with the rotation or how they are defined.
Looking at the local stiffness matrix K (before rotation), a few components are not equal, but they should be.
It should be that
K(5,3) = - K(6,2)
but they are not equal
(this applies to other components as well, which follows from symmetry)

   2100000         0         0         0         0         0 ...
         0     20360         0         0         0      9566 ...
         0         0     20360         0    -10180         0 ...
         0         0         0      1346         0         0 ...
         0         0    -10180         0      6841         0 ...
         0      9566         0         0         0      6841 ...
...

So, the source for those terms
K_ij = B^T_ik D_km B_mj
which reduces (diagonal D matrix in this case)
K_62 = B^T_6k D_kk B_k2 = B_k6 D_kk B_k2
K_53 = B^T_5k D_kk B_k3 = B_k5 D_kk B_k3
and the sum taking into account all the zero components we get the further simplfications:
K_62 = B_26 D_22 B_22 + B_66 D_66 B_62
K_53 = B_35 D_33 B_33 + B_55 D_55 B_53
where D_22 = D_33 and D_55 = D_66 (which I can see is fullfilled.
should be satisfied (and it is not).
Continue to dig and I end up with the core difference:
B_55 * B_53 == B_66 * B_62
These terms should be equal, but now, they differ in sign.
Note: There is other components as well due to symmetries

Changing the signs of
B_6,6 and B_6,12
results in a stiffness matrix that is identical to the beam2d case in both the x-z and x-y planes.
Assuming that someone checked the beam2d case carefully, I'm fairly certain this is the correct fix for this problem.

(the N matrix should probably change sign in the same positions, but I haven't checked closely).

Preferably, I'd like to find a resource which derives the B matrix used in this element. All the texts I find jump directly to a hardcoded K-matrix, which is no good (I'm quite frustrated that so many textbooks do this)

Hi, I was trying to have a look at the beam3d problem, but I am not sure if the problem is still there or has been resolved already.
To test it, I have created a model of simply supported beam with square cross section. As the bending around y and z principal axes is independent, I should get the same displacements and rotations. The height/length ration is 0.1/5 =  0.2, so definitely a thin beam. Beam is loaded by a constant load.
The deflection in the middle should be 5/384*(f*l^4/EI) = 2,929804692, rotation at the supports should be fl^3/(24/*EI)=1,875075, so at x=0 the rot_y = -1,875075, rot_z=1,875075 and at the end (x=5) rot_y=1,875075 and rot_z=-1,875075.
The local end forces reported by individual elements are also ok.
And this is what I get (input file attached) from the repo (oofem.org) version.

Thank you all for the support.
Please have a look to the attached files: in a cantilever beam, when I put:

beamShearCoeff 0.8333333333333

the element still misses the equilibrium. If you try with

beamShearCoeff 0.8333333333333

it's ok except but the moment at the free end is not zero (see here).
Is the Timoshenk's contribution really supported?

thanks

And does that output file really match the input file in the zip?
The output file shows deflection in both y and z.
Your input file is only loaded in z-direction, and when i run it, i only get deflection in z (like expected).

I also looked at the beamshearcoeff (I'm not to familiar with these things), but I cannot see it being used, well, anywhere.
Instead, the only place I find it is here

    //shearCoeff = this->give(CS_BeamShearCoeff);
    shearAreay = this->give(CS_ShearAreaY, gp);
    shearAreaz = this->give(CS_ShearAreaZ, gp);

which makes it clear that BeamShearCoeff has been replaced, and these 2 areas are used instead (but someone forgot to update the documentation).

Ah, my mistake (i thought i tested the two versions and found the same output, but I must have been looking at the wrong files.
The shear coefficient is indeed used to compute the shear areas.
As in  "shear_area_z = shear_area_y = beamshearcoeff * area;"
I'm not at all familiar with these shear coefficients, so I'm just going with what I see in the code (though, it seems I easily miss things).

Thanks Mikael.
Ok, but it is definitely complicated, at least in my opinion.the post processor needs to know the shape functions...
I say again that we must find a way to account for the external load even if it is applied via sets, the output is pretty confusing now.

I made a simple example to illustrate what i meant about other element types.

Lets have a bar prescribed at the left hand side, and free at the right hand side, and we will discretize this to a single element.

|[                 ] Free end


Distributed load along x-direction:
|[-> -> -> -> -> ->]


Projected load onto FE-approximation:
|[                 ]---->


Integration points along the bar element:
|[  x     x     x  ]

All three integration points will contain the same stress, as this problem is from the point of view of the element, is indistinguishable from that of a concentrated nodal force.
ted nodal load.

The same thing happens with all element types (2D/3D), and the only difference is that we don't have an analytical value to compare against, as those are not statically determined.
The only option then is to refine the problem (higher order or finer mesh).

Perhaps this is a language issue (I'm not a native english speaker myself), but I mean statically determined as opposed to statically indeterminate. In the 1D specimens, there is only one point where the reaction forces apply, as opposed to the distributed of line forces where you cut out a piece of a plate.
So all I tried to say was that this type of post-processing is only possible for 1D specimens.

I only wanted to point out that there is nothing special going on with the beam here; The integration point values always reflect a load that is projected on the shape functions. The shear force matches the beam kinematics and the shear coefficient the user supplies. So there is no lies in the output. If there is a fine enough mesh, the integration-point moments will match the analytical solution. But for a coarse mesh, it'll be way off, just like with the bar.

When it comes to local end forces, yes, the output misses loads applies by sets, and this is certainly not intentional. It's a side-effect of the way the standard output writes data, that there is no (nice) way to get to the loads from sets. This might seem like I'm being a tad stubborn, but it will complicate quite a few plans I had.
I'd be much better if we got rid of the "local end forces" from the output, and made a specialized beam export tool.
(This reminds me much of how ANSYS workbench deals with beams, where you add the "beam tool" to the results)