% MotionGenesis file: SphericalPendulumConstraintsSoft.txt % Problem: Pendulum on a string (wrecking ball or tether ball) % Copyright (c) 2009 Motion Genesis LLC. All rights reserved. %-------------------------------------------------------------------- NewtonianFrame N Particle Q %-------------------------------------------------------------------- Variable x'', y'', z'' % Nx>, Ny>, Nz> measures of Q's position from No. Constant g = 9.8 m/s^2 % Earth's gravitational acceleration. Constant k = 100 N/m % Linear spring constant (when string is flexible). Constant zeta = 0 noUnits % For linear damping constant (when string is flexible). Constant bAir = 0 N*s/m % Aerodynamic damping constant. Constant Ln = 0.3 m % Natural length of spring. Constant epsilonL = 1.0E-17 m % Avoid divide-by-zero error when L is near zero. Q.SetMass( m = 1 kg ) %-------------------------------------------------------------------- % Position, velocity, acceleration Q.Translate( No, x*Nx> + y*Ny> + z*Nz> ) %-------------------------------------------------------------------- % L, L', L'' are useful for both constraints and forces. % L = Q.GetDistance(No) + epsilonL. Specified L+, L', L'' SolveQuadraticPositiveRootDt( L^2 - Q.GetDistanceSquared(No), L ) L += epsilonL %-------------------------------------------------------------------- % String's spring tension on Q (spring can stretch, but cannot compress). % In other words, the spring is active when SpringStretch > 0. Specified SpringStretch' % Spring stretch and its time-derivative. SpringStretch = L - Ln isSpringActive = ceil( sign(SpringStretch)/2 ) TensionSpring = isSpringActive * k * SpringStretch QsUnitVectorFromNo> = Q.GetPosition(No) / L FSpring> = -TensionSpring * QsUnitVectorFromNo> %-------------------------------------------------------------------- % String's damper force on Q, active when SpringStretch' > 0. % Damper is active only when spring is active and stretching. % Note: To facilitate picking a damping constant, use damping ratio % with zeta=0 for no damping, zeta=1 for critical damping. b = 2 * zeta * sqrt(m*k) SpringStretch' = Dt( SpringStretch ) isElongating = ceil( sign(SpringStretch')/2 ) isDamperActive = isSpringActive * isElongating TensionDamper = isDamperActive * b * SpringStretch' FDamper> = -TensionDamper * QsUnitVectorFromNo> %-------------------------------------------------------------------- % Air's damping force on Q is -bAir * magV * (v / magV). FAir> = -bAir * Q.GetVelocity(N) %-------------------------------------------------------------------- % Rewrite net force on Q in terms of Fx, Fy, Fz. FNet> = -m*g*Ny> + FAir> + FSpring> + FDamper> Fx = Dot( FNet>, Nx> ); Fy = Dot( FNet>, Ny> ); Fz = Dot( FNet>, Nz> ) FNet> := Fx*Nx> + Fy*Ny> + Fz*Nz> Q.AddForce( FNet> ) %-------------------------------------------------------------------- % Equations of motion. Zero[1] = Dot( Nx>, Q.GetDynamics() ) Zero[2] = Dot( Ny>, Q.GetDynamics() ) Zero[3] = Dot( Nz>, Q.GetDynamics() ) %-------------------------------------------------------------------- % Track kinetic energy, gravitational energy, spring energy. KE = Q.GetKineticEnergy() PEGravity = System.GetForceGravityPotentialEnergy( -g*Ny>, No ) PESpring = 1/2 * k * SpringStretch^2 * isSpringActive %-------------------------------------------------------------------- % Calculate power and work associated with damping forces. PowerDamper = Dot( FDamper>, Q.GetVelocity(N) ) PowerAir = Dot( FAir>, Q.GetVelocity(N) ) Variable WorkDamper' = PowerDamper Variable WorkAir' = PowerAir %-------------------------------------------------------------------- % Kinetic energy + potential energy - work done should remain constant. EnergySum = KE + PEGravity + PESpring - WorkDamper - WorkAir %-------------------------------------------------------------------- % Calculate other output quantities. Tension = TensionSpring + TensionDamper q = GetAngleBetweenUnitVectors( -Ny>, QsUnitVectorFromNo> ) %-------------------------------------------------------------------- % Integration parameters and initial values for work. Input tFinal = 16 sec, tStep = 0.005 sec, absError := 1.0E-07, relError = 1.0E-07 Input WorkAir = 0 Joules, WorkDamper = 0 Joules %-------------------------------------------------------------------- % Default (zero) initial values for variables (override below). Input x = 0 m/s, y = 0 m/s, z = 0 m/s Input x' = 0 m/s, y' = 0 m/s, z' = 0 m/s %-------------------------------------------------------------------- % Initial values for variables correspond to swinging spring. % Input x := Input(Ln) * sinDegrees(1) m, y := -Input(Ln) * cosDegrees(1) m % OutputPlot t sec, L m % OutputPlot t sec, q degrees %-------------------------------------------------------------------- % Alternate input values to swing from -90 deg to 150 deg and free-fall. Input tFinal := 2.0 sec, tStep := 0.001 sec, k := 10000 N/m, zeta := 1.0 noUnits Input x := -Input(Ln) m, y' := -2.35 m/s OutputPlot x m, y m OutputPlot t sec, Tension Newton %-------------------------------------------------------------------- % Alternate input values to spiral like a tether-ball. % Input tFinal := 16.0 sec, tStep := 0.001 sec, k := 5.0E6 N/m, zeta := 1.0 noUnits, bAir := 0.5 N*sec/m % Input x := -Input(Ln) m, z' := -70 m/s % OutputPlot x m, z m % OutputPlot t sec, KE Joules, WorkAir Joules, EnergySum Joules %-------------------------------------------------------------------- % List output quantities and solve ODEs. Output t sec, q deg, L m, SpringStretch m, L' m/s, L'' m/s, x m, y m, z m, Tension Newton, & KE Joules, PEGravity Joules, PESpring Joules, WorkDamper Joules, WorkAir Joules, EnergySum Joules ODE( Zero, x'', y'', z'' ) SphericalPendulumConstraintsSoft %-------------------------------------------------------------------- Save SphericalPendulumConstraintsSoft.html Quit