Correct Displacements, but incorrect moments using the DSM for Structure

[Tygra_1983]

Dear all,

I have the following structure with its loading that I am trying to solve



I am using the stiffness method, and my MATLAB code is below:

clear, clc, close all

% Structural properties
E = 2.1E+08;
I = 1022000/100^4;
A = 620/100^2;
EI = E*I
EA = E*A
alpha = atand(30/40)
L1 = 10;
L2 = 40;
L3 = 50;
L4 = 30;
N = 12;

% Local stiffness matrix for member 1
K1 = [EA/L1 0 0 -EA/L1 0 0;
      0 12*EI/L1^3 6*EI/L1^2 0 -12*EI/L1^3 6*EI/L1^2;
      0 6*EI/L1^2 4*EI/L1 0 -6*EI/L1^2 2*EI/L1;
      -EA/L1 0 0 EA/L1 0 0;
      0 -12*EI/L1^3 -6*EI/L1^2 0 12*EI/L1^3 -6*EI/L1^2;
      0 6*EI/L1^2 2*EI/L1 0 -6*EI/L1^2 4*EI/L1;]

% Local stiffness matrix for member 2
K2 = [EA/L2 0 0 -EA/L2 0 0;
      0 12*EI/L2^3 6*EI/L2^2 0 -12*EI/L2^3 6*EI/L2^2;
      0 6*EI/L2^2 4*EI/L2 0 -6*EI/L2^2 2*EI/L2;
      -EA/L2 0 0 EA/L2 0 0;
      0 -12*EI/L2^3 -6*EI/L2^2 0 12*EI/L2^3 -6*EI/L2^2;
      0 6*EI/L2^2 2*EI/L2 0 -6*EI/L2^2 4*EI/L2;]

% Local stiffness matrix for member 3
K3 = [EA/L3 0 0 -EA/L3 0 0;
      0 12*EI/L3^3 6*EI/L3^2 0 -12*EI/L3^3 6*EI/L3^2;
      0 6*EI/L3^2 4*EI/L3 0 -6*EI/L3^2 2*EI/L3;
      -EA/L3 0 0 EA/L3 0 0;
      0 -12*EI/L3^3 -6*EI/L3^2 0 12*EI/L3^3 -6*EI/L3^2;
      0 6*EI/L3^2 2*EI/L3 0 -6*EI/L3^2 4*EI/L3;]

% Local stiffness matrix for member 4
K4 = [EA/L4 0 0 -EA/L4 0 0;
      0 12*EI/L4^3 6*EI/L4^2 0 -12*EI/L4^3 6*EI/L4^2;
      0 6*EI/L4^2 4*EI/L4 0 -6*EI/L4^2 2*EI/L4;
      -EA/L4 0 0 EA/L4 0 0;
      0 -12*EI/L4^3 -6*EI/L4^2 0 12*EI/L4^3 -6*EI/L4^2;
      0 6*EI/L4^2 2*EI/L4 0 -6*EI/L4^2 4*EI/L4;]


% Transformation matrix for member 1
T1 = zeros(6,N);
T1(1,1) = 1;
T1(2,2) = 1;
T1(3,3) = 1;
T1(4,4) = 1;
T1(5,5) = 1;
T1(6,6) = 1;

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

% Transformation matrix for member 3
T3 = zeros(6,N);
T3(1,4) = cosd(alpha);
T3(1,5) = sind(alpha);
T3(3,6) = 0;
T3(4,10) = cosd(alpha);
T3(4,11) = sind(alpha);
T3(6,12) = 0

% Transformation matrix for member 4
T4 = zeros(6,N);
T4(2,7) = -1;
T4(1,8) = 1;
T4(3,9) = 1;
T4(5,10) = -1;
T4(4,11) = 1;
T4(6,12) = 1

% Mapping from local to global using transformation matrix
Km1 = T1'*K1*T1;
Km2 = T2'*K2*T2;
Km3 = T3'*K3*T3;
Km4 = T4'*K4*T4;

% Sum up all element global stiffness matrices
Km = Km1 + Km2 + Km3 + Km4

% Application of boundary conditions
Km(:,[1,2,8]) = []
Km([1,2,8],:) = []

% Strucrtural load vector
F = zeros(N,1);
F(2) = 150;
F(3) = 250;
F(5) = 150 + 600
F(6) = -250 + 4000
F(8) = 600;
F(9) = -4000

% Application of boundary conditions
F([1,2,8]) = []

% Solving for displacements in metres
U = inv(Km)*F

% Entire displacement vector
Ux = [0;0;U(1:5);0;U(6:9)]

Ux =

         0
         0
    0.0600
    0.0000
    0.5479
    0.0449
   -0.0005
         0
   -0.0377
    0.4122
   -0.0003
   -0.0018



% Find internal forces and moments

Ax = T2*Ux % For member 2
A = K2*Ax

A =

   1.0e+04 *

    0.0171            % Axial Force
    0.0278            % Shear
    1.0000            % Bending Moment
   -0.0171            % Axial Force
   -0.0278            % Shear
    0.1137            % Bending Moment

Strangely, I am getting the correct displacements at the nodes, but when I try to retrieve the internal forces and moments I am getting incorrect answers. However, I am getting the correct axial forces, but the shear and moment are incorrect. Please observe vector A at the end of the code.

Please see, the bending moments for the structure below:




This is quite strange. Does anyone know why the moments I am getting are incorrect?

Many thanks,

Tygra

 

[rb1957]

I'm not familiar with this s/ware but ...

1) can you explain your load vector ? your loading looks to be distributed loads, but your load vector looks to be point loads.
2) do you mean you "weld" element 2 and 4 together at N2 ?
3) what constraint at N3 ? I'd've thought X and Y, but it looks like the Y load is being reacted by element 4 ?? ... is N3 a constraint ? or are elements 3 and 4 a "frame" on the beam ?/ (if so, why any load ??)
4) what do your boundary condition statements mean ? where is N8 ?
5) is "E" correct ? You're using N m units, looks like I and A are corrected from cm, E = 2.1E8 Pa (I don't work in metric)

 

[Tygra_1983]

Hi rb1957,

1) So, these are the fixed end moments. the UDL is 30 kN/m. Thus, for member 1 that is 10m long, the force at the node is 30*10/2 = 150 kN, and the moment is 30*10^2/12 250 kNm. For member 2 that is 40m long, the force at the node is 30*40/2 = 600 kN, and the moment is 30*40^/12 = 4000 kNm. Therefore, at node 1 we have 150 kN and 250 kNm. At node 2 we have (150kN + 600kN) and (-250kNm + 4000kNm). At node 3 we have 600N and -4000kNm.

2) Yes, element 2 and element 4 are welded and fixed. Thus, there is a bending moment.

3) Element 3 is pinned at both ends. Thus, it only transmits axial force. Element 1, 2 and 4 are fixed.

4) For each node there are 3 degrees of freedom: vertically, horizontally and rotationally. There are 4 nodes in total, so 4*3 = 12 degrees of freedom for the entire structure. At node 3 this means that the horizontal is DOF 7, Vertical is DOF 8 and Rotation is DOF 9. At this node there is a roller support so the vertical is restrained. Hence, why N8 is a boundary condition and removed from the global stiffness matrix.

5) I am using kN m units, and E is 2.1E+08 kPa.

I hope this helps.

 

[rb1957]

1) wow, what a complicated way to apply load ! "At node 3 we have 600N and -4000kNm." ... you mean at node 5 ? So F(2) is FY at N1, and F(3) is MZ at N1 ? N1 (the LH node on element 1) has +ve MZ and N2 (the RH node on element 2) has +ve MZ (as you've written it) ... but no F(8) in FY at N2 and F(9) is MZ at N2 (F(5) and F(6) have two loads therefore N5 ... yes ?).

4) ok, so there are only three reactions, ok, fine so ...
a) N3 is just floating, element 3 offers some limited support to the beam.
b) If "8" is FY at Node 2 then you've built the displacement matrix as FX, FY, MZ at N1, N5, N2, N3 (or N1, N3, N2, N5 ... but the former looks more natural ?).

You can see you need to explain the model in a lot more detail so we can understand it.

 

[Tygra_1983]

Thanks for you reply, rb1957. After much staring at the screen and wondering why (lol) I have realised that you have to add on the fixed end forces/moments to the vector A in the code.

Thus, for member 2 we get

Member2_forces = [
A(1)
A(2) + 600
A(3) + 4000
A(4)
A(5) + 600
A(6) - 4000]

Member2_forces =   1.0e+03 *  

-0.1712    0.3216   -6.0000    0.1712    0.8784   -5.1374

These are the correct values!

 

[rb1957]

glad to hear you solved your problem ! after the "palmslap" there's a feeling of solving something !

 

[Tygra_1983]

Thanks rb1957.

 

[rb1957]

if you want some more "pain", try solving with the unit force method, let the load in element 3 be the redundant force (remove element 3 and the resulting structure is statically determinate). As I think of it, the solution maybe be trickier than normal ... the typical redundancy is in the Y direction, not inclined ... I wonder how this'll change the solution ? N5 would deflect in the Y direction, so maybe it is how does this 5Y deflection load up members 3 and 4 (as opposed to a deflection along element 3)?