Why Are My Bending Moment Signs Incorrect Using the Direct Stiffness Method?

[Tygra]

Homework Statement

Why am I getting incorrect signage on my structure

Relevant Equations

Kglobal = T^T*KLocal*T

Hi there,

I was wondering if someone might know why I am getting incorrect signs for my structure that I am working on using the Direct Stiffness Method? I am following the procedure that I was taught when I was at University and I can't completely remember everything. I am designing a seven-storey structure that looks like this:

The magnitude of the forces I am getting is accurate, but the signs are incorrect.

Here are the bending moment on the structure:

You might not be able to see, but for the columns, you get a positive bending moment at the bottom of each storey and a negative bending moment at the top of the storey - this pattern is the same for each storey.

In my code I am getting positive bending moments at the top and the bottom, but like I said the magnitude is quite accurate.

Rather than post the entirity of my code, lets focus on a single column.

The local stiffness matrix for a column has the form:

And I set up the transformation matrix as follows:


Code:

T =

T =

  6x6 table

              U1    V1    theta1    U2    V2    theta2
              __    __    ______    __    __    ______

    U1         0    1       0        0    0       0
    V1        -1    0       0        0    0       0
    theta1     0    0       1        0    0       0
    U2         0    0       0        0    1       0
    V2         0    0       0       -1    0       0
    theta2     0    0       0        0     0      1

To compute the global stiffness matrix you use the equation:

This will give you the global stiffness matrix! You can then proceed to find the displacements and rotations. Then you can find the internal forces on the structure.

I guess the thing to look at is the transformation matrix. Does it seem correct to you. The local stiffness matrix is definitely correct.

Many thanks in advance!

 

[Lnewqban]

What makes the structure deform?
How rigidly are the links joined?

 

[Arjan82]

Please specify the problem a bit more clearly. What are the loads you are applying?

Can you say why you think it is only a minus sign? That is not clear to me at all, also because I don't know which loads you are applying.

The colors are not that clear. Red positive, blue negative? X-axis positive to the right, y-axis positive upwards? That means a positive moment is counter clockwise?

I think the rotation matrix rotates the beam clockwise by 90 degrees, is that what you expected?

 

[Tygra]

Sorry for the ambiguity, guys.

I am applying horizontal loads at the joints from left to right. The top load is 3.3 kN, the rest are 6.6 kN. The frame is a rigid frame with rigid joints, the bottom columns are fully fixed.

Purple is positive and blue is negative for the bending moments. Yes right is positive for the x axis and upward positive for the y axis.

The columns and beams bend in double curvature due to the shear. Thus, I should be getting a positive and negative bending moment at the column ends (that is what the software is giving me). However, I am getting positive bending moments at the top and bottom ends of the columns (using the stiffness method in MATLAB.

 

[Arjan82]

Ok, so the results of the bending moments are not your Matlab results, but that of some presumably correct software? If you don't get the negative part in the beam it seems more sensible to me that it is either a boundary condition problem (of the individual beams) or a stiffness matrix problem, because that means the force balance is not correct.

Maybe try a much simpler problem first? Just two vertical beans with one horizontal one, and just one load at one node?

 

[Tygra]

Hi Arjan, so I have done as you suggested a more simple structure. This time a 2 storey frame. I am still getting the same problem; accurate magnitudes, but incorrect signs. As this is a much smaller structure I will paste my code. There is a table at the end with the internal forces.

clear, clc, close all

% Structural Propeties
Ic = 22760/100^4;
E = 2.1E+08;
Ac = 170.2/100^2;
h = 4;
L = 6;

% Local stiffness matrix for columns
Kc = [E*Ac/h 0 0 -E*Ac/h 0 0;
      0 12*E*Ic/h^3 6*E*Ic/h^2 0 -12*E*Ic/h^3 6*E*Ic/h^2;
      0 6*E*Ic/h^2 4*E*Ic/h 0 -6*E*Ic/h^2 2*E*Ic/h;
      -E*Ac/h 0 0 E*Ac/h 0 0;
      0 -12*E*Ic/h^3 -6*E*Ic/h^2 0 12*E*Ic/h^3 -6*E*Ic/h^2;
      0 6*E*Ic/h^2 2*E*Ic/h 0 -6*E*Ic/h^2 4*E*Ic/h]

% Local Stiffness Matric for beams
Kb = [E*Ac/L 0 0 -E*Ac/L 0 0;
      0 12*E*Ic/L^3 6*E*Ic/L^2 0 -12*E*Ic/L^3 6*E*Ic/L^2;
      0 6*E*Ic/L^2 4*E*Ic/L 0 -6*E*Ic/L^2 2*E*Ic/L;
      -E*Ac/L 0 0 E*Ac/L 0 0;
      0 -12*E*Ic/L^3 -6*E*Ic/L^2 0 12*E*Ic/L^3 -6*E*Ic/L^2;
      0 6*E*Ic/L^2 2*E*Ic/L 0 -6*E*Ic/L^2 4*E*Ic/L]

% number of nodes
n = 18


% Transformation matrix for Member 1
T1 = zeros(6,n)
T1(2,1) = -1;
T1(1,2) = 1;
T1(3,3) = 1;
T1(5,4) = -1;
T1(4,5) = 1;
T1(6,6) = 1;

% Transformation matrix for Member 2
T2 = zeros(6,n)
T2(2,4) = -1;
T2(1,5) = 1;
T2(3,6) = 1;
T2(5,7) = -1;
T2(4,8) = 1;
T2(6,9) = 1;

% Transformation matrix for Member 3
T3 = zeros(6,n)
T3(2,10) = -1;
T3(1,11) = 1;
T3(3,12) = 1;
T3(5,13) = -1;
T3(4,14) = 1;
T3(6,15) = 1;

% Transformation matrix for Member 4
T4 = zeros(6,n)
T4(2,13) = -1;
T4(1,14) = 1;
T4(3,15) = 1;
T4(5,16) = -1;
T4(4,17) = 1;
T4(6,18) = 1;

% Transformation matrix for Member 5
T5 = zeros(6,n)
T5(1,4) = 1;
T5(2,5) = 1;
T5(3,6) = 1;
T5(4,13) = 1;
T5(5,14) = 1;
T5(6,15) = 1;

% Transformation matrix for Member 6
T6 = zeros(6,n)
T6(1,7) = 1;
T6(2,8) = 1;
T6(3,9) = 1;
T6(4,16) = 1;
T6(5,17) = 1;
T6(6,18) = 1;

% Combine transformation matrices into one vector
T = cat(3,T1,T2,T3,T4,T5,T6);

% Assemble global matrix
for i = 1:6
    Km(:,:,i) = T(:,:,i)'*Kc*T(:,:,i)
    if i > 4
        Km(:,:,i) = T(:,:,i)'*Kb*T(:,:,i)
    end
end

% sum up global matrices
Kt = Km(:,:,1)
for i = 1:5
    Kt = Kt + Km(:,:,i+1);
end

% Apply boundary conditions
Kt(:,[1,2,3,10,11,12]) = [];
Kt([1,2,3,10,11,12],:) = [];

% Assign 50 kN to each storey at their respective nodes
F = zeros(n,1);
F([4,7]) = 50;

% Apply the boundary conditions for the force vector
F([1,2,3,10,11,12]) = [];

% Solve for displacewments
U = inv(Kt)*F

% Re-construct a column vector containing all displacements
Ux = [0;0;0;U(1:6);0;0;0;U(7:12)]

% Compute internal forces for ecah member
for i = 1:6
    Fx(:,i) = T(:,:,i)*Ux
    Vx(:,i) = Kc*Fx(:,i)
    if i > 4
       Vx(:,i) = Kb*Fx(:,i)
    end
end

% set up table to display internal forces
rows = {'Axial 1','Shear 1','Moment 1','Axial 2','Shear 2','Moment 2'}
T = array2table(Vx,...
    'VariableNames',{'Member 1','Member 2','Member 3','Member 4','Member 5','Member 6'},'RowNames',rows)


T =

  6x6 table

                Member 1    Member 2    Member 3    Member 4    Member 5    Member 6
                ________    ________    ________    ________    ________    ________

    Axial 1     -57.396     -20.657      57.396      20.657      24.838       25.04
    Shear 1      50.121       24.96      49.879       25.04      -36.74     -20.657
    Moment 1      128.1       37.89      127.52       38.17     -110.28     -61.949
    Axial 2      57.396      20.657     -57.396     -20.657     -24.838      -25.04
    Shear 2     -50.121      -24.96     -49.879      -25.04       36.74      20.657
    Moment 2     72.388      61.949      71.991      61.991     -110.16     -61.991

Here are the results from the software:

I don't know why this is happening!

 

[Arjan82]

I will try to have a look tomorrow

 

[Tygra]

Okay, thank you.

 

[Arjan82]

So... that took me a bit longer... But I hadn't forgotten you :) I was confused as well, but I think I might know what is going on.

FIrst off, I think your code is entirely correct. The rotational matrix rotates a beam counter clockwise (contrary to what I said, because the transformation matrix is defined to transform global coordinates to local ones, not the other way around... Annoyingly I always make that mistake, even when I anticipate making that mistake I make that mistake....)
So, to be absolutely sure we are talking about the same things:
Element 1 is the lower left column. Its first node (with dofs no 1,2,3) is attached to the ground, its second node is the top one (with dofs no 4,5,6) connecting to the rest of the structure. What you've calculated in the table at the end of the code is the forces on this element in the local coordinates (i.e. positive axial is upwards, positive shear is to the left, positive rotation is ccw)

On node 2 of element 1 there is a horizontal force to the right (i.e. negative shear from the element's perspective). So on node 1 there needs to be a positive moment (i.e. ccw) to compensate for this force. Which is indeed what you calculate.

If the column was not attached to anything else, the column would just bent and the direction normal to the column would be somewhat downward. Due to the rest of the structure (horizontal beam 5, and vertical beam 3) this is however not possible and you need to apply another positive moment to bent the beam 'back' to vertical.

In other words, you indeed get two positive moments on the nodes. You can apply the same reasoning to all other columns and beams.

So, what about the software? Well, maybe it defines this as postive:

Note that the positive direction of the moments now switch sign at each of the nodes, and with this definition you would indeed have a negative moment at node 2 of element 1. This would be consistent with the results of your software.

 

[Tygra]

Hi Arjan, sorry for my very late reply as well. I actually posted the same question on the Engineering Tips Forum. There was no activity on it so I posted here, and eventually someone responded to me over on Engineering Tips.

You are absolutely right, and the answer to the whole dilema is a very simple one: The direct stiffness results are the applied moments to the end of the member the internal moment, the plotted moment diagram from tekla (The software), comes from the application of statics/mechanics along the length of the member (Celt83, Eng Tips).

Therefore, there is one more step after retrieving the results from the Direct Stiffness Method, where you simply apply the mechanics/statics to the joint translations/rotations found, and then you will get the internal member data - this is what the software is displaying.

Furthermore, just draw a free-body diagram of one of the members with the results obtained from the DSM and find the bending moment diagram - it will give you a positive and negative bending moments at each end.

Again, sorry that I took so long to post and thank you for all you help, Arjan!