A programming tutorial and a peek inside the programmer's mind

This blog will show how I create a program. I tried to put down the thought process as I wrote the code.

Scilab is a free scientific programming language / IDE / toolbox similar to Matlab. But free!

This Scilab program will calculate the shear forces distributed to a fastener pattern by a load in the same plan as the fasteners, but whose location can be different from the location of the centroid of the fastener pattern. Direct shear and axial loads are distributed only as a function of the fastener material and diameter, not the fastener location.

The moment resulting from the eccentricity of the applied force and any applied couple moments are also distributed to the fasteners as shear forces.

 The theory behind the solution of such problems can be found in many books or articles.

To start, I pasted the 'Moments of Inertia' code into Scilab because it has a similar structure to what I am trying to accomplish. I went through the variables initialization and input section of the older program and changed names and erased or new variables as needed to fit the current problem.

Sometimes I realize I need more variables and I have to revisit this section of the code. Rather than comment a lot, I make the variable names as suggestive as possible. I believe it helps people who read this code not to have to refer to comments all the time.

Here is the code in its first iteration. Let’s see if this works:

____________________________________________________________________________________________

fasteners=evstr(input("Number of fasteners?"));

x=zeros(10);

y=zeros(10);

diam=zeros(10);

ShearU=zeros(10);

area=zeros(10);

DiamSq=zeros(10);

DiamSqX=zeros(10);

DiamSqY=zeros(10);

rSquare=zeros(10);

SumShearUrSquare=0;

ShearUr=zeros(10);

ShearUr2=zeros(10);

ShearX=zeros(10);

ShearY=zeros(10);

ShearMomentX=zeros(10);

ShearMomentY=zeros(10);

SumDiamSq=0;

SumDiamSqX=0;

SumDiamSqY=0;

Xcentroid=0;

Ycentroid=0;

SumSU=0;

Mreaction=0;

Preaction=0;

Vreacgtion=0;

P=evstr(input("Axial Force Component?"));

V=evstr(input("Shear Force Component?"));

forcex=(input("X location of force?"));

forcey=(input("Y location of force?"));

for i=1:fasteners

printf("for fastener %d \n",i);

x(i)=evstr(input("Input x coordinate of fastener center"));

y(i)=evstr(input("Input y coordinate of fastener center"));

diam(i)=evstr(input("Input diameter of fastener"));

ShearU(i)=evstr(input("Input Ultimate Shear Strength of fastener"));

 

The variables needed to calculate the centroid will be computed during this cycle as well.

 

DiamSq(i)=diam(i)^2;

DiamSqX(i)=DiamSq(i)*x(i);

DiamSqY(i)=DiamSq(i)*y(i);

 

As these variables are calculated, they are added to summing variables

 

SumDiamSq=SumDiamSq+DiamSq(i);

SumDiamSqX=SumDiamSqX+DiamSqX(i);

SumDiamSqY=SumDiamSqY+DiamSqY(i);

 

Reviewing the formulas, we see that the sum of the ultimate shear strengths is needed to distribute the shear and axial loads to the fasteners.

 

SumSU=SumSU+ShearU(i);

end

 

At the end of this first for cycle we have enough information to calculate the position of the centroid.

 

Xcentroid=SumDiamSqX/SumDiamSq;

Ycentroid=SumDiamSqY/SumDiamSq;

 

Another 'for' cycle is needed to calculate auxiliary variables used to distribute the shear forces to each fastener. r is the position vector from the centroid to each fastener.

for i=1:fasteners

rSquare(i)=(x(i)-Xcentroid)^2+(y(i)-Ycentroid)^2;

 

the products of the ultimate shear strength and r squared will also have to be summed. Quite a lot of boring coding here :)

 

SumShearUrSquare=SumShearUrSquare+rSquare(i)*ShearU(i);

end

 

The eccentricity of the applied load will generate a reaction moment at the centroid of the fastener pattern. Also, the force components will reverse direction (these are forces and moments transmitted to fasteners).

A little thinking here: if the problem were solved by hand, we could assume positive directions for forces and moments. They would be drawn as arrows on the FBD, and if calculations gave negative results the direction of the arrows on the FBD would be reversed. The program cannot see the FBD, thus we need to keep the sign until the very end. Also, any moment computation will have to be done through the cross product of the force and moment vectors.

First we need to determine the shear and axial forces at the centroid of the fastener pattern, and the moment generated by the eccentricity of the applied load. For the moment, the right hand rule convention is used, going into the 'board' (or clockwise).

Not being able to infer right away which sign to place on the variables V and P to get the correct direction for the reaction moment, I used an example problem in the book and set the signs to give me the same result.

 

Mreaction=V*(forcex-Xcentroid)+P*(forcey-Ycentroid);

Preaction=-P;

Vreaction=-V;

 

Now we can finally distribute these loads to the fasteners! Note that the sign of the forces will have to be AGAIN reversed to show forces acting ON the fasteners. Rather unnecessary, but hey, I'm just following the book. A better approach might be to have a clear idea of the solution BEFORE you start coding. I prefer to get directly to the program, rather than draw flowcharts or write pseudocode.

 

for i=1:fasteners

ShearX(i)=-Preaction*ShearU(i)/SumSU;

ShearY(i)=-Vreaction*ShearU(i)/SumSU;

 

For the moment created fastener shear forces, we have to again make sure that the sign we assign to each variable will generate a forces pointing in the correct direction according to the direction of the moment at the centroid of the pattern. We can do that by visually inspecting an example problem already worked by hand.

 

ShearMomentX(i)=-Mreaction*ShearU(i)*(y(i)-Ycentroid)/SumShearUrSquare;

ShearMomentY(i)=Mreaction*ShearU(i)*(x(i)-Xcentroid)/SumShearUrSquare;

 

The X and Y components from the direct forces and the moment are summed and the results displayed

 

ShearX(i)=ShearX(i)+ShearMomentX(i);

ShearY(i)=ShearY(i)+ShearMomentY(i);

printf('Fastener %d incurs a horizontal shear force of %f and a vertical shear force of %f\n',i,ShearX(i),ShearY(i));

end

 ___________________________________________________________________________________________

To correct for syntax errors...I just run the program and make changes until I get no errors. This is the first stage of debugging. The second stage will be to compare the results with 'the book' and make more changes as needed. The code might change, sometimes significantly, so I will paste it here again. Scilab is case sensitive....this might take some time.

Once all the syntax errors have been corrected and we're actually getting results, it's time for phase two of the debugging process: making sure the values we get are correct.

For that we will use a solved example problem and see if we get the same results as the solved example used for reference.

And here is the final result, which yields correct values for the distributed shear forces on each fastener. This is final iteration will go into the program library at http://www.solvengineer.com/codes.html

Scroll down to the end to see how the code was debugged.

___________________________________________________________________________________________ 

fasteners=evstr(input("Number of fasteners?"));

x=zeros(10);

y=zeros(10);

diam=zeros(10);

ShearU=zeros(10);

area=zeros(10);

DiamSq=zeros(10);

DiamSqX=zeros(10);

DiamSqY=zeros(10);

rSquare=zeros(10);

SumShearUrSquare=0;

ShearUr=zeros(10);

ShearUr2=zeros(10);

ShearX=zeros(10);

ShearY=zeros(10);

ShearMomentX=zeros(10);

ShearMomentY=zeros(10);

SumDiamSq=0;

SumDiamSqX=0;

SumDiamSqY=0;

Xcentroid=0;

Ycentroid=0;

SumSU=0;

Mreaction=0;

Preaction=0;

Vreacgtion=0;

AppliedMoment=0;

P=evstr(input("Axial Force Component?"));

V=evstr(input("Shear Force Component?"));

forcex=evstr(input("X location of force?"));

forcey=evstr(input("Y location of force?"));

AppliedMoment=evstr(input("Applied Moment? (Right Hand Rule sign convention going in is positive)"));

for i=1:fasteners

printf("for fastener %d \n",i);

x(i)=evstr(input("Input x coordinate of fastener center"));

y(i)=evstr(input("Input y coordinate of fastener center"));

diam(i)=evstr(input("Input diameter of fastener"));

ShearU(i)=evstr(input("Input Ultimate Shear Strength of fastener"));

DiamSq(i)=diam(i)^2;

DiamSqX(i)=DiamSq(i)*x(i);

DiamSqY(i)=DiamSq(i)*y(i);

SumDiamSq=SumDiamSq+DiamSq(i);

SumDiamSqX=SumDiamSqX+DiamSqX(i);

SumDiamSqY=SumDiamSqY+DiamSqY(i);

SumSU=SumSU+ShearU(i);

end

Xcentroid=SumDiamSqX/SumDiamSq;

Ycentroid=SumDiamSqY/SumDiamSq;

for i=1:fasteners

rSquare(i)=(x(i)-Xcentroid)^2+(y(i)-Ycentroid)^2;

SumShearUrSquare=SumShearUrSquare+rSquare(i)*ShearU(i);

end

Mreaction=-P*(forcey-Ycentroid)+V*(forcex-Xcentroid)-AppliedMoment;

Preaction=-P;

Vreaction=-V;

for i=1:fasteners

ShearX(i)=-Preaction*ShearU(i)/SumSU;

ShearY(i)=-Vreaction*ShearU(i)/SumSU;

ShearMomentX(i)=-Mreaction*ShearU(i)*(y(i)-Ycentroid)/SumShearUrSquare;

ShearMomentY(i)=Mreaction*ShearU(i)*(x(i)-Xcentroid)/SumShearUrSquare;

printf("Fastener %d incurs an axial force (horizontal load component) of %f\n",i,ShearX(i));

ShearX(i)=ShearX(i)+ShearMomentX(i);

printf("Fastener %d incurs a shear force (vertical load component) of %f\n",i,ShearY(i));

ShearY(i)=ShearY(i)+ShearMomentY(i);

printf("Fastener %d incurs a horizontal load from the centroid moment of %f and a vertical load from the centroid moment of %f\n",i,ShearMomentX(i),ShearMomentY(i));

printf('Fastener %d incurs a horizontal shear force of %f and a vertical shear force of %f\n',i,ShearX(i),ShearY(i));

end

printf("Fastener %d incurs an axial force (horizontal load component) of %f\n",i,ShearX(i));

ShearX(i)=ShearX(i)+ShearMomentX(i);

printf("Fastener %d incurs a shear force (vertical load component) of %f\n",i,ShearY(i));

ShearY(i)=ShearY(i)+ShearMomentY(i);

printf("Fastener %d incurs a horizontal load from the centroid moment of %f and a vertical load from the centroid moment of %f\n",i,ShearMomentX(i),ShearMomentY(i));

printf('Fastener %d incurs a horizontal shear force of %f and a vertical shear force of %f\n',i,ShearX(i),ShearY(i));

end

____________________________________________________________________________________________ 

After all the syntax error have been corrected, we run the program again and compare the fastener shear forces with the results from a solved problem. Unfortunately, our results are WRONG! Not a big surprise, that's what happens 99 times out of 100. Let's start fixing.

Running the program again we discover that the fastener forces from the horizontal and vertical loads (P and V) are correct. This means the moment at the centroid has not been distributed correctly to the fastener pattern.

First we need to verify if the sign of the moment-distributed loads is correct. Comparing the first fastener with the known results shows both the horizontal and vertical components are reversed.

This might give us a starting point in debugging the coded moment equations! However, the magnitudes are also incorrect so most likely there are multiple errors in the program.

An easy way to determine if at least the computed moment at the centroid is correct is to run the program again and find out the value of the moment from the command line (no need to load the code with displaying even more variables who are not needed)

The reaction moment - remember, this has the opposite sign as the moment generating the loads in the fasteners, is found to have the opposite sign as we established by the RHR, and is different in magnitude as well.

Mreaction=V*(forcex-Xcentroid)+P*(forcey-Ycentroid);

Looking at the formula, it seems X centroid and Y centroid have been miscalculated somehow. We can find their values easily from the command line. But it turns out Xcentroid and Ycentroid are correct..... 

So much for an easy fix :P We need to go over the way Mreaction was calculated again. We try to match the formula shown in the hand calculation of  M reaction, keeping in mind that the distances between the applied loads and the centroid might be negative numbers. If we get the formula right for this instance, it will work for any case. (hopefully).No, for real. You can verify yourselves.

Looks like we made a mistake the first time. Let's try this: Mreaction=-P*(forcey-Ycentroid)+V*(forcex-Xcentroid);.

Running the program again, we finally get the right results. Further testing is necessary to ensure this program works for all cases. In order to make the program apply to more real-life situations, we are going to add a couple moment to the applied forces.

The changes in the code should be miminal: AppliedMoment=evstr(input("Applied Moment? (Right Hand Rule sign convention going in is positive)"));

When we calculate the reaction moment the AppliedMoment will be added with a negative sign. Mreaction=-P*(forcey-Ycentroid)+V*(forcex-Xcentroid)-AppliedMoment; based on the RHR convention established.

Testing on yet another solved problem yields the correct answers!

Note: this program assumed that the direct axial and shear loads are distributed equally among similar fasteners (the only differences resulting from the different ultimate strengths), regardless of their position. In more realistic analysis scenarios, outer fasteners on a splice joint would be loaded more than the inner fasteners.

I hope this blog will help the scientific programmer understand some of the basics of solving real-world problems using the greatest automation tool ever created, the computer.

Thank you

Mihai Pruna