Torsion Waves—Solved Numerically

Basic System Properties

The actual system from Pasco has 75 rods, so set numberOfRods=75, but it can be any positive integer.

numberOfRods = 75;

κ = 1600.0 and I = 1.0 makes waves that go 40 rods per second in the long wavelength limit. This has not been calibrated against the actual system’s rods and connecting wires.

kappa = 1600.0;

moment = 1.0;

damping = 0.4;

Torque Calculation

To get the torque on a rod, we first need to capture the fact that each rod is connected by wires to two neighbors (except for the end ones, which we’ll get to in a minute).

neighbors[i_]:= {i -1, i + 1}

The angle of a rod relative to its neighbors has to be special-cased for the end rods. We do this by calculating the twist relative to the “far left neighbor at 0” and the “far right neighbor at 75” specially.

clampedLeft = 0.0; clampedRight = 0.0;

twist[rodAngle_,i_,neighbor_,t_]:=Which[neighbor==0, clampedLeft rodAngle[i][t], neighbor==numberOfRods+1,clampedRight rodAngle[i][t], True, rodAngle[i][t]-rodAngle[neighbor][t]]

If you want the left or the right end to be free, set clampedLeft or clampedRight to 0. If you want them clamped, you can set them to 1.

torque[rodAngle_, i_, t_]:=Sum[-kappa * twist[rodAngle, i, j, t], {j, neighbors[i]}]

Equations of Motion

rodEquations = Table[moment rodAngle[k]''[t]== torque[rodAngle, k, t]-damping *rodAngle[k]'[t], {k, 1, numberOfRods}];

Initial Conditions

An odd function centered at the center of the rods is our initial impulse (angular velocities):

impulseCenter = (numberOfRods + 1 )/2;
impulseWidth = 7.0;
impulseAmplitude = 1.0;
impulseFunction[k_]:=impulseAmplitude(k-impulseCenter) Exp[-(((k-impulseCenter))^2)/(impulseWidth^2)]

rodInitialVelocities = Table[rodAngle[k]'[0]==impulseFunction[k], {k, 1, numberOfRods}];

Our initial positions are zero:

rodInitialPositions = Table[rodAngle[k][0]==0, {k, 1, numberOfRods}];

Use NDSolve

Have Mathematica numerically solve the system for a duration of 15 seconds.

duration = 15.0;

allEquations=Join[rodEquations, rodInitialVelocities, rodInitialPositions];

allAngles =Table[rodAngle[k],{k, 1, numberOfRods}];

interpolatedSolutions = NDSolve[allEquations, allAngles, {t, 0, duration}];

Animate the Solution

Try watching the animation at 1/4 speed (so 15 seconds lasts one minute).

animationRate=1.0;

maxPlotAngle=Pi/6.0;

Animate[ListPlot[Evaluate[Table[Evaluate[{k,rodAngle[k][t]}],{k, 1, numberOfRods}]/.interpolatedSolutions],PlotRange->{-maxPlotAngle,maxPlotAngle}], {t, 0, duration}, AnimationRate->animationRate, AnimationRepetitions->1]

Notes

The dispersion relation for long wavelengths is (ω^2)=(κ/I)(((k s))^2), where s is the spacing between the rods, κ is the restoring torque for the wires connecting them, and I is the moment of inertia of a rod. However, this is an approximation and the actual dispersion relation is (ω^2)=(κ/I)(2-2cos(k s)). Therefore it is normal (not a numerical artifact) that the wave pulses change shape as the animation progresses. This can be seen in the actual Pasco system.

The dispersion is much more obvious if one turns the damping coefficient off or down.