Laboratory No. 2

System Dynamics
Carter Schultz in collaboration with Ethan Oberg

Task 1: Springs

Part 1: Determine the theoretical spring constant for the cantilever beam,

Our measurements of the cantilever beam were:

Blength = Quantity[4.0,"in"];

Bwidth = Quantity[.756, "in"];

Bthickness = Quantity[.167, "in"]; (*The Quantity command defines a value with units in Mathematica this can be used for automatic units*)

A quick Internet search yielded a young’s Modulus of 6061 Aluminum. http://asm.matweb.com/search/SpecificMaterial.asp?bassnum=MA6061t6

Eal= Quantity[10000,"ksi"];

The relationship between these quantities and the spring constant (k) is:

k == ((3 e i)/(L^3));

where, i is the area moment of inertia:

i == ((b (h^3))/12);

The spring constant can then be solved for as:

kCant = ((3*Eal*((Bwidth*(Bthickness^3))/12))/(Blength^3))

With correct significant figures this number is:

kCant = NumberForm[kCant,3](*This command displays the number to a defined precision *)

Part 2: Experimentally determine the spring constant

To experimentally determine the spring constant we devised an experiment with the following procedure. Weights were hung from the spring and the magnitude of the weight and deflection of the beam measured. The data produced from our experiment for the cantilever beam are:

Weight ={0,246.5,493.5,740.0,986.8,1234.0,986.8,740.0,493.5,246.5,0}*9.81/1000 ;(*Newtons*) 
Height = {10.29,10.21,10.092,9.943,9.875,9.754,9.868,9.993,10.104,10.216,10.243}; (*mm*)

Converting the height values to deflections:

Deflection = Height - 10.29

Arranging the values into a list of points

Plot1 = Table[{Weight[[i]],Deflection[[i]]},{i,11}];

Fitting a line to the data:

line1 = Fit[Plot1,{x},x]

So the equation (Displacement = -.0432 * Load) is the one term linear fit to the data.

Displaying both the data and line fit on the same axes

Show[ListPlot[Plot1,Mesh->All],Plot[{line1},{x,0,15}],AxesLabel->{"Load (Newtons)", "Displacement (mm)"},PlotLabel->"K Graph for Cantilever beam"]

Clearly the line fits the data very well so an experimentally determine k value is:

Kexperiement = -0.04317345643045399^-1

Kexperiement =Quantity[NumberForm[23.162379913010927 ,3],"N/mm"](*Newtons per mm*)

To compare this value to the predicted value I will convert it to the units of the prediction:

Kexperiement =NumberForm[N[UnitConvert[(QuantityUnit[23.2]),Quantity["kipf/in"]]],3]

Here I compare it to the predicted value by calculating the percentage error:

(.132-.138)/.138

Since there was less than 4% error between the prediction and experiment, I would call the exercise a great success.

Part 3: Energy analysis of Cantilever data

Numerically integrating the experimental data can yield the energy stored and released during the loading and unloading. This can be accomplished simply using Mathematica Integrate Function:

First I will split the list of points that represent the data into the loading and unloading sections:

loading =Plot1[[1;;6]]

unloading = Plot1[[6;;11]]

Then I will integrate both lists to get the energy of loading and unloading. Here I exploit Mathematica Interpolation function to preform a better than Trapezoidal approximation.

Integrate[Interpolation[loading,InterpolationOrder->3][x],{x,Max[loading[[All,1]]],Min[loading[[All,1]]]}]

Integrate[Interpolation[unloading,InterpolationOrder->3][x],{x,Max[unloading[[All,1]]],Min[unloading[[All,1]]]}]

So according to the data 3.14 mJ of energy were used to load the beam, and 3.05 mJ were used to unload it.

Part 4: Analysis of Complex Beam

The complex beam was testing using an identical procedure to the simple beam. The complex beam data was as follows:

Weight ={0, 246.5,493.5,740,986.8,1234,986.8,740,493.5,246.5,0}*9.81/1000 (*Newtons*);

Deflection = {12.985,12.606,12.222,11.849,11.468,11.092,11.462,11.835,12.214,12.601,12.986}-12.985; (*mm*)

To begin my analysis I will first preform a linear regression on the data and display the data an regression on the same plot:

Plot1 = Table[{Weight[[i]],Deflection[[i]]},{i,11}];

line1 = Fit[Plot1,{x},x]

So the equation (Weight = -.15707 * Deflection) best fits the data.

Show[ListPlot[Plot1,Mesh->All],Plot[{line1},{x,0,15}],AxesLabel->{"Load (Newtons)", "Displacement (mm)"},PlotLabel->"K Graph for Complex beam"]

The linear regression fits the data extremely well. Therefore, the coefficient from the regression can be used to express the k value of the spring.

-0.15707029705200887^-1

KvalueForComplexBeam = 6.37 (*N/mm*);

Interestingly even though it has extremely complex geometry the complex bean still has a perfectly linear response. While it would be extremely difficult to predict the k value of the complex beam, a simple experiment can very accurately predict the behavior of the system. Energy analysis of the complex beam is as follows:

First I split the data into two sets, one containing the loading points, and one containing the unloading points.

loading =Plot1[[1;;6]]

unloading = Plot1[[6;;11]]

Now I numerically integrate these values to determine the energy stored.

Integrate[Interpolation[loading,InterpolationOrder->3][x],{x,Max[loading[[All,1]]],Min[loading[[All,1]]]}]

Integrate[Interpolation[unloading,InterpolationOrder->3][x],{x,Max[unloading[[All,1]]],Min[unloading[[All,1]]]}]

According to the data 11.49 mJ of energy was used to load the beam, and 11.57 mJ recovered when unloading it. However, this is physically impossible and therefore the error bounds on our calculations must be greater than the differences in these two energies. It is likely that a single point being misread would be enough to cause this level of error, so human error is very likely.

Part 5: Analysis of Bellows Spring

The bellows spring data with a sealed vent was:

Weight = {0, 246.5,493.5,740,986.8,1234,986.8,740,493.5,246.5,0}*9.81/1000 (*Newtons*);

Deflection = {12.685,11.7,10.734,9.724,8.83,7.969,8.815,10.236,10.677,11.609,12.612} - 12.685 (*mm*);

Preforming a linear regression:

Plot1 = Table[{Weight[[i]],Deflection[[i]]},{i,11}];

line1 = Fit[Plot1,{x},x]

The produced linear equation is:

Show[ListPlot[Plot1,Mesh->All],Plot[{line1},{x,0,15}],AxesLabel->{"Load (Newtons)", "Displacement (mm)"},PlotLabel->"K Graph for Bellows Spring"]

Besides a single data point which clearly appears as an outlier, the bellows spring with a sealed vent almost perfectly fits the linear regression. No hysteresis appears in the data.

.392209^-1

K = 2.54966 N/mm for unvented bellows spring.

An energy analysis of the spring is as follows:

Separate the measured points into two data sets:

loading =Plot1[[1;;6]]

unloading = Plot1[[6;;11]]

Numerically integrate the two data sets to yield a resulting value for the energy:

Integrate[Interpolation[loading,InterpolationOrder->3][x],{x,Max[loading[[All,1]]],Min[loading[[All,1]]]}]

Integrate[Interpolation[unloading,InterpolationOrder->3][x],{x,Max[unloading[[All,1]]],Min[unloading[[All,1]]]}]

So 29.4 mJ were used to load the spring and 28.8 mJ were recovered. It appears there were some slight mechanical losses.

The Bellows Spring with an open vent was:

Deflection = {12.363,10.238,8.165,6.36,4.513,2.785,4.505,6.299,8.477,10.212,12.357}-12.363 (*mm*);

Weight = {0, 246.5,493.5,740,986.8,1234,986.8,740,493.5,246.5,0}*9.81/1000 (*Newtons*);

Preforming a linear regression:

Plot1 = Table[{Weight[[i]],Deflection[[i]]},{i,11}];

line1 = Fit[Plot1,{x},x]

Show[ListPlot[Plot1,Mesh->All],Plot[{line1},{x,0,15}],AxesLabel->{"Load (Newtons)", "Displacement (mm)"},PlotLabel->"K Graph for Bellows Spring(Vented)"]

Once again, the spring exhibited extremely linear behavior.

-.813526^-1

K = 1.22 N/mm for the vented bellows spring

loading =Plot1[[1;;6]]

unloading = Plot1[[6;;11]]

Integrate[Interpolation[loading,InterpolationOrder->3][x],{x,Max[loading[[All,1]]],Min[loading[[All,1]]]}]

Integrate[Interpolation[unloading,InterpolationOrder->3][x],{x,Max[unloading[[All,1]]],Min[unloading[[All,1]]]}]

60.56 mJ spend loading the beam,and 60.14 mJ recovered unloading it. I expected a greater energy loss for the vented bellows than the unvented, do to the air escaping from the vent experiencing viscous friction. However, the presence of this effect  does not appear within the data. Unvented the spring still exhibited linear behavior. However, the spring constant was significantly reduced.

Part 2: Dampers

To determine whether the supplied dash-pot was linear we devised an experiment. Our experiment involved moving the dash-pot to its upper end of travel, connecting a fixed weight to it and measuring the time it took the dash-pot to reach the other end of travel. This method was reproduced for a set of different weights. The data produced by the experiment is:

Weight= {22.72,44.92,67.21,89.39,112.09,134.61}*9.81/1000 (*Newtons*) ; Time = {33.5,18.6,13.9,10.5,8.6,7.5} (*Seconds*); Height= 3.33 (*in*);

Plotting the data yields:

Plot1 = Table[{Weight[[i]],Time[[i]]},{i,6}];

ListPlot[Plot1,Mesh->All, AxesLabel->{"Load (Newtons)", "Duration (s)"},PlotLabel->"Time for a given weight"]

This graph shows that there is a polynomial relationship between weight and time. However, this alone does not indicate whether or not the dash-pot is linear. Instead we need some information about the dynamics of the system.

When the system was observed there was no noticeable acceleration of the mass. We will assume that the amount of time the system took to accelerate was extremely short, and therefore the velocity profile can be considered constant. This will produce a slight overestimate for the magnitude of the damping, however it should not significantly effect results. Assuming near instant acceleration the steady state of the system is expressed as a force balance: weight == b* velocity -> weight == b*distance/time.

Solve[Weight=={b1,b2,b3,b4,b5,b6}*3.33/Time,{b1,b2,b3,b4,b5,b6}]

Plotting the b value over the weight yields:

bs ={b1,b2,b3,b4,b5,b6}/.{b1->2.2422183783783787`,b2->2.4613731891891897`,b3->2.752158675675676`,b4->2.7650501351351355`,b5->2.8398152972972976`,b6->2.974153378378379`}

Plot1 = Table[{Weight[[i]],bs[[i]]},{i,6}];

ListPlot[Plot1,Mesh->All,AxesLabel->{"Load (Newtons)", "Damping (N/(in/s))"},PlotLabel->"B Value over Weight"]

Clearly the b value for this spring does not appear constant under our instant acceleration assumption. Either our assumption is invalid or the damper is non-linear.

A mean value for the damping coefficient is:

Mean[bs]

So b can be approximated as 2.67 N/(in/s)

Finally the energy dissipated by the damper needs to be calculated. This energy could be solved for by predicting the force the damper would exert based on its b value and integrating this over the distance. However, I will exploit the conservation of energy:

The change in energy for the system is the change in potential energy from the initial state to the final state:

U==m*g*h;

U== weight *h;

U= Quantity[Weight,"N"] *Quantity[Height,"in"]

Converting to correct units and sig figs:

NumberForm[UnitConvert[U,Quantity[1,"mJ"]],3]

So the energy dissipation of the system ranges from 18.9 mJ to 112 mJ.