## Direct Stiffness Method Help

## Direct Stiffness Method Help

(OP)

Hi all,

I am trying to solve for a gable frame using the Direct Stiffness Method. The frame is obviously constructed of two rafters and two columns. I am quite confused: I am getting the correct results for the columns, but not for the rafters. I am comparing my results to those in ETABS for confirmation of the correct results. I am using my beloved Octave, so I need someone who uses either Octave or MATLAB.

The frame consists of two 8 metre columns, and two 15.083 metre rafters at an angle of 6 degrees to the horizontal. The UDL on the rafters is 13.14 kN/m

Here is the code from Octave:

The force vector gives the axial force, shear and bending moment for a single member over two joints. There are three degress of freedom at each joint.

This is the results for the column:

These results are correct!!

These are the results for the rafter:

These results are incorrect!

I have examined everything and cannot see what is the issue. I know this question is quite time comsuming, so I deeply thank anyone that attempts to help!

I am trying to solve for a gable frame using the Direct Stiffness Method. The frame is obviously constructed of two rafters and two columns. I am quite confused: I am getting the correct results for the columns, but not for the rafters. I am comparing my results to those in ETABS for confirmation of the correct results. I am using my beloved Octave, so I need someone who uses either Octave or MATLAB.

The frame consists of two 8 metre columns, and two 15.083 metre rafters at an angle of 6 degrees to the horizontal. The UDL on the rafters is 13.14 kN/m

Here is the code from Octave:

#### CODE -->

clear, clc, close all% Structural dataE = 2.1E+08; Ic = 8.74999E-04;Second Moment areas of ColumnsAc = 0.01442;Cross-sectional areas of columnsLc = 8; n = 15; alpha = 6Slope of roof in degreesmatrix = zeros(6,n); Lr = 15/(cosd(alpha)) Ir = 2.94331E-04Second Moment areas of raftersAr = 0.00856;Cross-sectional areas of rafters% Local stiffness matrix for columnsKc = [1 0 0 -1 0 0 0 12 6*Lc 0 -12 6*Lc 0 6*Lc 4*Lc^2 0 -6*Lc 2*Lc^2 -1 0 0 1 0 0 0 -12 -6*Lc 0 12 -6*Lc 0 6*Lc 2*Lc^2 0 -6*Lc 4*Lc^2]% Local stiffness matrix for raftersKr = [1 0 0 -1 0 0 0 12 6*Lr 0 -12 6*Lr 0 6*Lr 4*Lr^2 0 -6*Lr 2*Lr^2 -1 0 0 1 0 0 0 -12 -6*Lr 0 12 -6*Lr 0 6*Lr 2*Lr^2 0 -6*Lr 4*Lr^2]% For loop to complete the local stiffness matricesfor i = 1:6 if i == 1 || i == 4 Kc(i,:) = Kc(i,:)*(E*Ac)/Lc; Kr(i,:) = Kr(i,:)*(E*Ar)/Lr; else Kc(i,:) = Kc(i,:)*(E*Ic)/Lc^3; Kr(i,:) = Kr(i,:)*(E*Ir)/Lr^3; end end% Transformation matricesT1 = matrix;For Left ColumnT1(2,1) = -1; T1(1,2) = 1; T1(3,3) = 1; T1(5,4) = -1; T1(4,5) = 1; T1(6,6) = 1 T4 = matrix;For right columnT4(2,13) = -1; T4(1,14) = 1; T4(3,15) = 1; T4(5,10) = -1; T4(4,11) = 1; T4(6,12) = 1; T2 = matrix; T2(1,4) = cosd(alpha);For left rafterT2(1,5) = sind(alpha); T2(2,4) = -sind(alpha); T2(2,5) = cosd(alpha); T2(3,6) = 1; T2(4,7) = cosd(alpha); T2(4,8) = sind(alpha); T2(5,7) = -sind(alpha); T2(5,8) = cosd(alpha); T2(6,9) = 1; T3 = matrix; T3(1,7) = cosd(-alpha);For right RafterT3(1,8) = sind(-alpha); T3(2,7) = -sind(-alpha); T3(2,8) = cosd(-alpha); T3(3,9) = 1; T3(4,10) = cosd(-alpha); T3(4,11) = sind(-alpha); T3(5,10) = -sind(-alpha); T3(5,11) = cosd(-alpha); T3(6,12) = 1;% Assemblykm1 = T1'*Kc*T1; km2 = T2'*Kr*T2; km3 = T3'*Kr*T3; km4 = T4'*Kc*T4; Ks = km1 + km2 + km3 + km4% Application of boundary conditionsKs(:,[1,2,13,14]) =[]; Ks([1,2,13,14],:) =[];% Applying loadingw = 13.14; V = w*Lr/2; M = w*Lr^2/12; F = zeros(n,1); F(4) = -V*sind(alpha); F(5) = V*cosd(alpha); F(6) = M; F(8) = 2*V*cosd(alpha); F(9) = 0; F(10) = V*sind(alpha); F(11) = V*cosd(alpha) F(12) = -M; F([1,2,13,14]) = []% SolvingU = inv(Ks)*F Ux = [0 0 U(1:10)' 0 0 U(11)']; fx = T2*Ux' force = Kr*fx

The force vector gives the axial force, shear and bending moment for a single member over two joints. There are three degress of freedom at each joint.

This is the results for the column:

#### CODE -->

Ux = [0 0 U(1:10)' 0 0 U(11)']; fx = T1*Ux' force = Kc*fx force = -197.1000 108.3198 0 197.1000 -108.3198 866.5586

These results are correct!!

These are the results for the rafter:

#### CODE -->

Ux = [0 0 U(1:10)' 0 0 U(11)']; fx = T2*Ux' force = Kr*fx force = -128.329 -85.605 -617.462 128.329 85.605 -673.685

These results are incorrect!

I have examined everything and cannot see what is the issue. I know this question is quite time comsuming, so I deeply thank anyone that attempts to help!

## RE: Direct Stiffness Method Help

T1 (left column)

## CODE -->

T4 (Right column)

## CODE -->

T2 (left rafter)

## CODE -->

T3 (right rafter)

## CODE -->

## RE: Direct Stiffness Method Help

I'd try the following, because it seems like there's a good chance your rotation transformation might be the problem.

Single horizontal beam with axial force and end moments.

Single vertical beam with axial force and end moments.

Then step up to a two member frame with one vertical and one horizontal member. (Fix the ends of one member to cut down the DOFs.) This one is to test whether your member DOFs area mapping correctly to your structure DOFs.

## RE: Direct Stiffness Method Help

First is consistent units. There is a lot going on the DSM and plenty of places to mess up units. Ensure that you have properly converted all units as needed before the matrix solution begins.

Next question, when you say the results are wrong are we talking order of magnitude wrong or just a little bit off?

SAP2000 has the ability to include shear deflections which you are not doing here. So if you want a reasonable comparison ensure those shear deflections are neglected in SAP.

Finally you have used a zero stiffness assignment for your boundary conditions, ensure that in the solution phase you are properly partitioning the global stiffness matrix into restrained and unrestrained DOFs. Its easy for something to be lost in translation here.

## RE: Direct Stiffness Method Help

the load on a column should be 13.14*15*cosd(6) = 196 kN ...

interestingly, 197 kN is 15*13.14 ... the inclined load on the roof ?

the six results for the column are three forces at either end, yes? why if force3 unbalanced

nah, more like 2 inplane forces and the inplane moment ... 108*8 = 864 ... ok, but that is a large lateral load on the column ...

surely the applied lateral load is 13.14*15*sind(6) ? wouldn't this be divided between the two ends of the roof ? (so the peak sees equal and opposite loads ?)

"Hoffen wir mal, dass alles gut geht !"

General Paulus, Nov 1942, outside Stalingrad after the launch of Operation Uranus.

## RE: Direct Stiffness Method Help

## RE: Direct Stiffness Method Help

I tried a single cantilever beam at an angle and got correct results. At first I thought the transformation matrix for the rotation would be the culprit, but now I don't think it is.

We are talking order of magnitude. For the rafter, they results are way off.

My results:

## CODE -->

Results in ETABS:

Could you elborate a bit more please, Celt83? Do you mean at node 3 (central node) F(9) of the force vector?

## RE: Direct Stiffness Method Help

K U = F + Ffixed

To recover F:

F = K U - Ffixed

Your Etabs result seems to indicate the axial force (F1) in the rafter changes over itโs length so either self weight is being considered or you have applied the load on the horizontal projection instead of perpendicular to the rafter local axis per your sketch.

## RE: Direct Stiffness Method Help

## RE: Direct Stiffness Method Help

## RE: Direct Stiffness Method Help

I have another question if thats okay.

In ETABS, I am trying to replicate an Eaves Haunch. The way I have done this was to combine to seperate cross sections. One that tapers and one that is straight. Obviously, this creates a node joint between the members. So, my question is, will this combined member act as one single member, or is the joint (even though it is a rigid joint) not allow for the two seperate members to act as one single member? Because in the real world they weld a cut member to the bottom of the rafter to create the eaves haunch - not two seperate members.

Here is what it looks like:

Thanks again!

## RE: Direct Stiffness Method Help

Where you need to be careful is in design because the member length will typically default to be node to node so this will impact unbraced lengths for column/beam designs. Typically software like Etabs, Stadd, etc will allow you to manually define Lb for each axis of the member to override the computed node-to-node distance.

Here are my results for a simple fixed-free beam with a point load at the end (N6 and N8):

Here are my inputs:

Units are: inches for length and Kips for load

Sign Convention:

Positive loads are in the direction of the load axis

Positive moments are counter-clockwise

x+ to the right โ

y+ is upward ๐

## CODE --> python

## RE: Direct Stiffness Method Help

"Hoffen wir mal, dass alles gut geht !"

General Paulus, Nov 1942, outside Stalingrad after the launch of Operation Uranus.

## RE: Direct Stiffness Method Help

Try an L-shaped frame with vertical column and horizontal beam as the next step. Fix the frame at both supports to cut the DOfs down to three. That problem is small enough -- unlike your gable frame -- so that you could practically develop the manual solution and know exactly where your m file is going wrong.

## RE: Direct Stiffness Method Help

I am now considering second order effects. However, I'm not sure whether to calculate the factor by which you amplify the forces by, either by the wind load on the structure, the notional horizontal force, or the combination of both.

## RE: Direct Stiffness Method Help

Do you have a manual approximation that you can use to check that your program is working correctly? The B1, B2 method is great for this. If not, then you know by this point in the thread what I'd recommend. LOL!!

## RE: Direct Stiffness Method Help

Hi 271828,

I am using the Direct Stiffness Method in Octave. Changing the subject slightly, I have modelled the eaves haunch for the gable by cutting it into five equal segments, so in reality it would have a staggered appearance (not a tapered one). The haunch is about 3m long, so, each segement is 0.6m. Obviously, the second moment of area decreases as you move towards the sharp end of the haunch, because the cross-section height decreases. In terms of results I am getting 4.51mm deflection from a force in ETABS, and 4.01mm using the same force using DSM in Octave. Do you think my method and results are good approximations?

Many thanks.

## RE: Direct Stiffness Method Help

## RE: Direct Stiffness Method Help

Most stuff I work on 0.5mm is an insignificant amount. For your specific application you would need to judge that for yourself.

Seems like you modeled the tapered haunch the right way, I do the same approach.

## RE: Direct Stiffness Method Help

## RE: Direct Stiffness Method Help

If you want to model the stiffness in the panel zone, then I'd recommend trying either a krawinkler or scissors approach.

Simple running the member properties all the way over to the node usually does an ok job of modeling flexibility within the panel zone.