MGBeamOnTwoCablesDynamicsViaFBD.html  (MotionGenesis input/output).
   (1) %    File: MGBeamOnTwoCablesDynamicsViaFBD.txt
   (2) % Problem: Dynamics of a beam supported by two cables.
   (3) %--------------------------------------------------------------------
   (4) NewtonianFrame N         % Ceiling with Nx> horizontally right, Ny> vertically down
   (5) RigidBody      B         % Beam with Bx> pointing from Bo to Bc and Bz> = Nz>.
   (6) Point          Nc( N )   % Point of N attached to cable C.
   (7) Point          Bc( B )   % Point of B attached to cable C.
   (8) %--------------------------------------------------------------------
   (9) Constant  LN = 6 m       % Distance between No and NC.
   (10) Constant  LB = 4 m       % Distance between Bo and BC.
   (11) Constant  LA = 2.7 m     % Length of cable A (connects Bo to No).
   (12) Constant  LC = 3.7 m     % Length of cable C (connects BC to NC).
   (13) Constant  g = 9.8 m/s^2  % Earth's gravitational acceleration.
   (14) Constant  b = 230 N*s/m  % Aerodynamic damping constant.
   (15) B.SetMass( m = 120 kg )
   (16) B.SetInertia( Bcm, Izz = m*LB^2/12, 0, Izz )
-> (17) Izz = 0.08333333*m*LB^2

   (18) %--------------------------------------------------------------------
   (19) Variable  x'', y''       % Nx> and Ny> measures of Bo's position vector from No.
   (20) Variable  q''            % Angle from Nx> to Bx> with positive sense +Nz>.
   (21) Variable  TA, TC         % Tensions in cables A and C.
   (22) %--------------------------------------------------------------------
   (23) %   Rotational kinematics relating  Bx>, By>, Bz>  to  Nx>, Ny>, Nz>.
   (24) B.RotateZ( N,  q )
-> (25) B_N = [cos(q), sin(q), 0;  -sin(q), cos(q), 0;  0, 0, 1]
-> (26) w_B_N> = q'*Bz>
-> (27) alf_B_N> = q''*Bz>

   (28) %--------------------------------------------------------------------
   (29) %   Translational kinematics.
   (30) Bo.Translate( No,  x*Nx> + y*Ny> )
-> (31) p_No_Bo> = x*Nx> + y*Ny>
-> (32) v_Bo_N> = x'*Nx> + y'*Ny>
-> (33) a_Bo_N> = x''*Nx> + y''*Ny>

   (34) Bcm.Translate( Bo,  0.5*LB*Bx> )
-> (35) p_Bo_Bcm> = 0.5*LB*Bx>
-> (36) v_Bcm_N> = 0.5*LB*q'*By> + x'*Nx> + y'*Ny>
-> (37) a_Bcm_N> = -0.5*LB*q'^2*Bx> + 0.5*LB*q''*By> + x''*Nx> + y''*Ny>

   (38) BC.SetPosition( Bo,  LB*Bx> )
-> (39) p_Bo_Bc> = LB*Bx>

   (40) NC.SetPosition( No,  LN*Nx> )
-> (41) p_No_Nc> = LN*Nx>

   (42) %--------------------------------------------------------------------
   (43) %   Constraint relating cable lengths to magnitudes of position vectors.
   (44) LengthConstraint[1] = LA^2 - Bo.GetDistanceSquared( No )
-> (45) LengthConstraint[1] = LA^2 - x^2 - y^2

   (46) LengthConstraint[2] = LC^2 - Bc.GetDistanceSquared( Nc )
-> (47) LengthConstraint[2] = LC^2 + 2*LB*cos(q)*(LN-x) - LB^2 - y^2 - 2*LB*y*sin(q)
        - (LN-x)^2

   (48) %--------------------------------------------------------------------
   (49) %   Forces on beam.
   (50) unitVectorBoToNo> = No.GetPositionVector(Bo) / LA
-> (51) unitVectorBoToNo> = -x/LA*Nx> - y/LA*Ny>

   (52) unitVectorBcToNc> = NC.GetPositionVector(Bc) / LC
-> (53) unitVectorBcToNc> = -LB/LC*Bx> + (LN-x)/LC*Nx> - y/LC*Ny>

   (54) Bo.AddForce( TA*unitVectorBoToNo> )
-> (55) Force_Bo> = -TA*x/LA*Nx> - TA*y/LA*Ny>

   (56) BC.AddForce( TC*unitVectorBcToNc> )
-> (57) Force_Bc> = -LB*TC/LC*Bx> + TC*(LN-x)/LC*Nx> - TC*y/LC*Ny>

   (58) Bcm.AddForce( m*g*Ny> - b * Bcm.GetVelocity(N) )
-> (59) Force_Bcm> = -0.5*b*LB*q'*By> - b*x'*Nx> + (m*g-b*y')*Ny>

   (60) %--------------------------------------------------------------------
   (61) %   Dynamics for B via FBD (free-body-diagram).
   (62) Dynamics[1] = Dot( Nx>,  B.GetDynamics()   )  % F = ma in Nx> direction.
-> (63) Dynamics[1] = TA*x/LA + b*x' + 0.5*LB*cos(q)*(2*TC/LC-m*q'^2) + m*x''
        - TC*(LN-x)/LC - 0.5*LB*sin(q)*(b*q'+m*q'')

   (64) Dynamics[2] = Dot( Ny>,  B.GetDynamics()   )  % F = ma in Ny> direction.
-> (65) Dynamics[2] = TA*y/LA + TC*y/LC + b*y' + 0.5*LB*sin(q)*(2*TC/LC-m*q'^2)
        + m*y'' + 0.5*LB*cos(q)*(b*q'+m*q'') - m*g

   (66) Dynamics[3] = Dot( Nz>,  B.GetDynamics(Bcm) )  % M = I*alpha in Nz> direction.
-> (67) Dynamics[3] = 0.5*LB*(TA*(x*sin(q)-y*cos(q))/LA+TC*(y*cos(q)+sin(q)*(
        LN-x))/LC) + Izz*q''

   (68) %--------------------------------------------------------------------
   (69) %   Integration parameters and initial values.
   (70) Input  tFinal = 9 sec, tStep = 0.02 sec,  absError = 1.0E-07
   (71) Input  x = 2 m,  x' = 0 m/s
   (72) %--------------------------------------------------------------------
   (73) %   Use length constraints to set initial values for y and q (with guess).
   (74) SolveSetInput( LengthConstraint = 0,   y = 3 m,  q = 15 deg )

->   %  INPUT has been assigned as follows:
->   %   y                         1.813835714721706       m
->   %   q                         27.67399398437716       deg

   (75) %--------------------------------------------------------------------
   (76) %   Differentiate length constraints to set initial values for y' and q' (with guess).
   (77) SolveSetInput( Dt(LengthConstraint) = 0,  y' = 0 m/s,  q' = 0 m/s )

->   %  INPUT has been assigned as follows:
->   %   y'                        0                       m/s
->   %   q'                        0                       m/s

   (78) %--------------------------------------------------------------------
   (79) %   Quantities to be output by ODE command.
   (80) Output  t sec,  x m,  y m,  q degrees,  TA Newtons,  TC Newtons
   (81) %--------------------------------------------------------------------
   (82) %   Append 2nd-derivative of LengthConstraint to dynamics equations and integrate.
   (83) ODE( [Dynamics; DtDt(LengthConstraint)] = 0,  x'', y'', q'', TA, TC )  MGBeamOnTwoCablesDynamics

   (84) %--------------------------------------------------------------------
   (85) %   Optional: Alternate 2-step method for solving Dynamics with Constraints.
   (86) %   1. Solve Dynamics[] equations for  x'', y'', q''  in terms of TA, TC.
   (87) %   2. Substitute these into the DtDt(LengthConstraint) to solve for TA, TC.
   (88) Solve( Dynamics = 0,   x'', y'', q'' )
-> (89) x'' = 0.5*(2*TC*(LN-x)/LC+b*LB*sin(q)*q'-2*TA*x/LA-2*b*x'-LB*cos(q)*(2*
        TC/LC-m*q'^2))/m - 0.25*LB^2*sin(q)*(TA*(x*sin(q)-y*cos(q))/LA+TC*(y*
        cos(q)+sin(q)*(LN-x))/LC)/Izz

-> (90) y'' = g + 0.25*LB^2*cos(q)*(TA*(x*sin(q)-y*cos(q))/LA+TC*(y*cos(q)+sin(
        q)*(LN-x))/LC)/Izz - 0.5*(2*TA*y/LA+2*TC*y/LC+2*b*y'+b*LB*cos(q)*q'+LB*
        sin(q)*(2*TC/LC-m*q'^2))/m

-> (91) q'' = -0.5*LB*(TA*(x*sin(q)-y*cos(q))/LA+TC*(y*cos(q)+sin(q)*(LN-x))/LC)
        /Izz

   (92) Solve( DtDt(LengthConstraint) = 0,   TA, TC )
-> (93) TA = 2*LA*((4*(LB^2+y*(y+2*LB*sin(q))+(LN-x)*(LN-x-2*LB*cos(q)))/m+LB^2
        *(sin(q)^2*(LN-x)^2+y*cos(q)*(y*cos(q)+2*sin(q)*(LN-x)))/Izz)*(2*x'^2+2
        *y'^2+x*(LB*cos(q)*q'^2-b*(2*x'-LB*sin(q)*q')/m)+y*(2*g+LB*sin(q)*q'^2-
        b*(2*y'+LB*cos(q)*q')/m))+(x*(4*(LN-x-LB*cos(q))/m-LB^2*sin(q)^2*(LN-x)
        /Izz)-y*(4*(y+LB*sin(q))/m-LB^2*cos(q)*(y*cos(q)+sin(q)*(LN-2*x))/Izz))
        *(2*x'^2+2*y'^2+LB^2*q'^2+4*LB*cos(q)*q'*y'+LB*cos(q)*(LN-x)*q'^2+2*LB*
        sin(q)*(g-b*y'/m)+y*(2*g-LB*sin(q)*q'^2-b*(2*y'+LB*cos(q)*q')/m)-4*LB*
        sin(q)*q'*x'-b*(2*LB*cos(q)*x'-(LN-x)*(2*x'-LB*sin(q)*q'))/m))/((y^2*(4
        /m+LB^2*cos(q)^2/Izz)+x*(4*x/m+LB^2*sin(q)*(x*sin(q)-2*y*cos(q))/Izz))*
        (4*(LB^2+y*(y+2*LB*sin(q))+(LN-x)*(LN-x-2*LB*cos(q)))/m+LB^2*(sin(q)^2*
        (LN-x)^2+y*cos(q)*(y*cos(q)+2*sin(q)*(LN-x)))/Izz)+(4*(y^2+LB*x*cos(q)+
        LB*y*sin(q)-x*(LN-x))/m+LB^2*(x*sin(q)-y*cos(q))*(y*cos(q)+sin(q)*(LN-x))
        /Izz)*(x*(4*(LN-x-LB*cos(q))/m-LB^2*sin(q)^2*(LN-x)/Izz)-y*(4*(y+LB*sin
        (q))/m-LB^2*cos(q)*(y*cos(q)+sin(q)*(LN-2*x))/Izz)))

-> (94) TC = -2*LC*((4*(y^2+LB*x*cos(q)+LB*y*sin(q)-x*(LN-x))/m+LB^2*(x*sin(q)-
        y*cos(q))*(y*cos(q)+sin(q)*(LN-x))/Izz)*(2*x'^2+2*y'^2+x*(LB*cos(q)*q'^2
        -b*(2*x'-LB*sin(q)*q')/m)+y*(2*g+LB*sin(q)*q'^2-b*(2*y'+LB*cos(q)*q')/m))
        -(y^2*(4/m+LB^2*cos(q)^2/Izz)+x*(4*x/m+LB^2*sin(q)*(x*sin(q)-2*y*cos(q))
        /Izz))*(2*x'^2+2*y'^2+LB^2*q'^2+4*LB*cos(q)*q'*y'+LB*cos(q)*(LN-x)*q'^2
        +2*LB*sin(q)*(g-b*y'/m)+y*(2*g-LB*sin(q)*q'^2-b*(2*y'+LB*cos(q)*q')/m)-
        4*LB*sin(q)*q'*x'-b*(2*LB*cos(q)*x'-(LN-x)*(2*x'-LB*sin(q)*q'))/m))/((y^2
        *(4/m+LB^2*cos(q)^2/Izz)+x*(4*x/m+LB^2*sin(q)*(x*sin(q)-2*y*cos(q))/Izz))
        *(4*(LB^2+y*(y+2*LB*sin(q))+(LN-x)*(LN-x-2*LB*cos(q)))/m+LB^2*(sin(q)^2
        *(LN-x)^2+y*cos(q)*(y*cos(q)+2*sin(q)*(LN-x)))/Izz)+(4*(y^2+LB*x*cos(q)
        +LB*y*sin(q)-x*(LN-x))/m+LB^2*(x*sin(q)-y*cos(q))*(y*cos(q)+sin(q)*(LN-
        x))/Izz)*(x*(4*(LN-x-LB*cos(q))/m-LB^2*sin(q)^2*(LN-x)/Izz)-y*(4*(y+LB*
        sin(q))/m-LB^2*cos(q)*(y*cos(q)+sin(q)*(LN-2*x))/Izz)))

   (95) %--------------------------------------------------------------------
   (96) %   Save input and output responses.
Saved by Motion Genesis LLC.   Command names and syntax: Copyright (c) 2009-2019 Motion Genesis LLC. All rights reserved.