The Torsional Oscillator --  Part 1

Clear["Global`*"]; DateString[]

Static Torque Measurements

Mathematica can’t do LinearModelFit on data with units, so we will  just take care to explicitly convert everything to the base SI units:  meters, kilograms, and seconds.  Enter the data as {angle, mass} pairs.  Call the mass positive if it tends to decrease the angle.

mvsθraw = (({{1.400, 450}, {1.505, 400}, {1.645, 350}, {1.800, 300}, {1.980, 250}, {2.180, 200}, {2.380, 150}, {2.600, 100}, {2.800, 50}, {3.005, 0}, {3.205, -50}, {3.405, -100}, {3.620, -150}, {3.830, -200}, {4.025, -250}, {4.235, -300}, {4.460, -350}, {4.625, -400}, {4.820, -450}, {3.005, 0}}));

We want the data in the form {τ, θraw}, which is backwards from the order in which it was entered above.  We also want to convert from hanging mass (in grams) to torque (in N-m).  This is easily done as a table.  The factor of ‘2’ is because there are 2 hanging masses.

The sign we assign to torque depends on whether we wish to consider the torque we apply, or the restoring torque (opposite sign) supplied by the wire.  It d

g = 9.8; (* m/s^2 *)
Raxle = QuantityMagnitude[
UnitConvert[
Quantity[1.00, "Inches"]/2, "Meters"]
]

θvsτ = Table[{2 * (mvsθraw[[i, 2]]/1000) * g * Raxle, mvsθraw[[i, 1]]}, {i, 1, Length[mvsθraw]}];

The high torque (or large angle) points tend to deviate a little from the line, so exclude them.  The exact parameters don’t matter much, but here I only include torques of absolute value < 0.1 Nm (which corresond to angular displacements of less than about 1.5 radians).

sfit = LinearModelFit[Select[θvsτ, -0.1 < #[[1]] < 0.1 &], {1, τ}, τ]

Show[{ListPlot[θvsτ, PlotStyle -> Red], Plot[sfit[τ], {τ, Min[θvsτ[[All, 1]]], Max[θvsτ[[All, 1]]]}]},
Frame -> True, FrameLabel -> {"Torque (N m)", "Angle (radians)"}, LabelStyle -> Larger, ImageSize -> Scaled[0.8]]

sfit["BestFitParameters"]

sfit["ParameterConfidenceIntervalTable"]

Look at the table entries to figure out how to extract the values and uncertainties.  θraw is the first line, and 1/κ is the second line.  Use the Around[ ] function to handle the uncertainties.  Store the fit table in a variable to reduce tedious retyping.

sparams = sfit["ParameterConfidenceIntervalTableEntries"];
MatrixForm[sparams]

θraw = Around[sparams[[1, 1]], sparams[[1, 2]]]

κstatic = ((-1)/(Around[sparams[[2, 1]], sparams[[2, 2]]]))

Dynamic Torque Measurements

Value of ΔI for added masses

Theoretical Derivation (not needed for the report)

Let ρ = the mass density of brass.  The quadrants are of height h.  Calculate Io for a solid disk of radius Ro, and Ii for a solid disk of radius Ri.  Subtracting them gives ΔI.  If Mo is the mass of the hypothetical solid outer disk, and Mi is the mass of the hypothetical solid inner disk, then the “M” given in the manual is the mass of the annular quadrant, 
M = Mo - Mi.

Clear[Mo, Mi, Ro, Ri, ρ, h, Io, Ii, ΔI, M]

Mo = ρ π (Ro^2)h;Io = (1/2)Mo (Ro^2) ; (* Mass and inertia for outer solid disk *)

Mi = ρ π (Ri^2)h;Ii = (1/2)Mi (Ri^2) ; (* Mass and inertia for inner solid disk *)

ΔI = Io - Ii; Simplify[ΔI]

ρ = ρ /. First[Solve[M == Mo - Mi, ρ]]

ΔI = Simplify[ΔI]

Numerical Values

M = Quantity[Around[0.2145 , 0.0005], "Kilograms"];
Ro = (1/2)Quantity[3.72, "Inches"];Ri = (1/2)Quantity[1.72, "Inches"];

Ri = UnitConvert[Ri, "Meters"]

Ro = UnitConvert[Ro, "Meters"]

ΔI = (1/2)M (Ro^2 + Ri^2)

ΔImag = QuantityMagnitude[ΔI]["Value"]

Enter the raw data for angular frequency as a function of number of added quadrants.

ωvsn = {{0, 5.397}, {2, 4.776}, {4, 4.331}, {6, 3.99}, {8, 3.721}};

Tvsn = Table[{ωvsn[[i, 1]], 2 π/ωvsn[[i, 2]]}, {i, 1, Length[ωvsn]}]

tvsn[i0_, n_, κ_] := 2 π  Sqrt[((i0 + n ΔImag)/κ) ]

dataplot = ListPlot[Tvsn, LabelStyle -> Larger, AxesLabel -> {"n", "Period"}, PlotStyle -> Red];

Clear[κ, i0]
dfit = NonlinearModelFit[Tvsn, tvsn[i0, n, κ],  {{i0, 0.001}, {κ, 0.06}}, n]

Show[
{Plot[dfit[n], {n, 0, 8}], dataplot}, 
Frame -> True, FrameLabel -> {"Number of added Brass Quadrants", "Period (s)"}, LabelStyle-> Larger, ImageSize -> Scaled[0.8]
]

dfit["ParameterConfidenceIntervalTable"]

dparams = dfit["ParameterConfidenceIntervalTableEntries"];
MatrixForm[dparams]

κdynamic = Around[dparams[[2, 1]], dparams[[2, 2]]]

Materials calculation for a predicted value of κ

Take values from section 6.1 of the instruction manual, included with the writeup.

See Eq. 3 in the writeup.  The extra 2 at the end is because there are really 2 wires, one on top and one on bottom, holding the rotor in place.  The material is “piano wire”, which is ASTM A228.  The shear modulus was taken from www.matweb.com.  Values in the instruction manual are typically diameters, and are also typically given in inches.  Include the uncertainty in the wire radius, since it is raised to the 4th power, that can be significant.

Mathematica can track and convert units, which is useful for this calculation.

κmat = ((π G (r^4))/(2L))* 2;

G = Quantity[80 (10^9), "Newtons/Meter^2"]

r = UnitConvert[
(1/2)Around[Quantity[0.039, "Inches"], Quantity[0.0005, "Inches"]],
"Meters"]

L = Quantity[0.254, "Meters"];

κmat

Comparison of results

TableForm[{ {QuantityMagnitude[κmat]}, {κstatic}, {κdynamic}},
TableHeadings->{{"Materials", "Static", "Dynamic"}, {"κ (N m)"}}]

The static and dynamic values disagree by just slightly more than one uncertainty, while the value estimated from the materials parameters has enough uncertainty to include either of the other two values.

κstatic - κdynamic

QuantityMagnitude[κmat] - κstatic

QuantityMagnitude[κmat] - κstatic == 0