How do I solve a symbolic system for a certain vector?

Hello everyone!

I've run into the rigid-body dynamics equation below:

(1)

which, reduced for a constant ω (i.e., ), results in:

(2)

In theory, for any inertia tensor J, there is a certain angular velocity ω that spins the rigid-body along its major axis of inertia, thus the spin is stable. Therefore, in that case, and

I've been able to prove it in a simulation environment using Simulink, but I want to find the theoretical relation between ω and τ, i.e., an equation system that takes equation (2) and writes it in terms of:

(3)

I've tried to expand it manually but I believe it is easier to solve using MATLAB Symbolic Toolbox. So far, I've been ably to write it for the thrust τ, i.e., simply using equation (2) and the following code snippet:

syms(sym('w',[1 3],'real'));
w = [w1; w2; w3];

syms(sym('T',[1 3],'real'));
T = [T1; T2; T3];

syms(sym('J',[3 3],'real'));
J = [J1_1 J1_2 J1_3; J2_1 J2_2 J2_3; J3_1 J3_2 J3_3];

T = cross(w,J*w);

My question is, is there any way to redefine the result from T = cross(w,J*w); and rewrite the equation in termos of ?

Thank you for your time!

Best regards,

Vítor

[Walter Roberson]

Note simpler code would be

syms w [1 3] real
syms T [1 3] real
syms J [3 3] real

You can see that w, T, J are already assigned the way you would want

w

w = 

T

T = 

J

J = 

[Vítor Narciso]

Duly noted, thanks! Any input on my issue?

[Walter Roberson]

r and D are not defined, I don't know what size they are

[Vítor Narciso]

Oopsy, that's a formula for J = R*D*R', but it does not matter for this context! Already corrected it!

[Sam Chak]

Hi @Vítor Narciso

To simplify the matter, if the inertia tensor in the body-fized frame is transformed to a diagonal matrix in the principal axes frame, called the Principal Moments of Inertia

then becomes

,

which implies that you want to apply the torques according to the formulas to make .

Thus, you want to solve these 3 equations simultaneously to obtain

[Vítor Narciso]

Yes, that's correct! In fact, I have solved by hand those equations for a diagonal inertia tensor (aka the Principal Moments of Inertia, as you stated). I wanted to use the full inertia tensor to solve it symbolically as, in my problem, I am dealing with Products of Inertia and it proves to be non-trivial to solve by hand.

Answer 1

[Sam Chak]

Hi @Vítor Narciso

It's an interesting problem as mentioned by @Paul. The following does not solve your problem symbolically in MATLAB, but an attempt to investigate your rigid body problem in two scenarios by the numerical approach.

To begin, run the following code accordingly:

tspan  = [0 10];
init   = [1; 2; 3];
opt = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);
[t, y] = ode45(@odefcn1, tspan, init, opt);      % Scenario 1
% [t, y] = ode45(@odefcn2, tspan, init, opt);    % Scenario 2
plot(t, y, 'linewidth', 1.5)
ylim([0 4]);

Scenario 1:

If the torque vector is designed as , then . Let's see what happens, but I guess you have already simulated this event in Simulink.

function dydt = odefcn1(t, y)

dydt = zeros(3,1);

J1 = 0.005; % principal moment of inertia about x-axis

J2 = 0.012; % principal moment of inertia about y-axis

J3 = 0.013; % principal moment of inertia about z-axis

dydt(1) = 0;

dydt(2) = 0;

dydt(3) = 0;

end

The angular velocities are stable (not diverging) but also not 'attractive' (not converging to the equilibrium point at ). It means that the rigid body has a stable rotation but tumbling. Can watch a video here:

https://www.youtube.com/watch?v=fZYfX9VvN8M

However, if the initial angular velocity is large, then this is generally undesirable, because the rigid body is at the risk of spinning and breaking up structurally weak components in the real world due to the centrifugal force.

More importantly, it is obeserved that when is applied to rigid body, the rotational motions are successfully decoupled and they become independent.

Scenario 2:

Thus, if the torque vector is designed as with , then and the angular velocities are expected to exponentially decay to zero from their non-zero initial angular velocities.

function dydt = odefcn2(t, y)

dydt = zeros(3,1);

J1 = 0.005; % principal moment of inertia about x-axis

J2 = 0.012; % principal moment of inertia about y-axis

J3 = 0.013; % principal moment of inertia about z-axis

k = 1;

tau1 = - (J2 - J3)*y(2)*y(3) - J1*k*y(1); % torque about x-axis

tau2 = - (J3 - J1)*y(3)*y(1) - J2*k*y(2); % torque about y-axis

tau3 = - (J1 - J2)*y(1)*y(2) - J3*k*y(3); % torque about z-axis

dydt(1) = ((J2 - J3)*y(2)*y(3) + tau1)/J1;

dydt(2) = ((J3 - J1)*y(3)*y(1) + tau2)/J2;

dydt(3) = ((J1 - J2)*y(1)*y(2) + tau3)/J3;

end

As theoretically predicted above, the response of the rigid body is exponentially stable and stop rotating at steady state.

Non-trivial Solutions:

This problem requires rigorous mathematical thinking to solve, at least to me. Need pencil and paper. Solving the 3 equations simultaneously yields two sets solutions (due to having 3 sets of quadratics). The solutions are valid for three non-zero torque components () and three different principal moments of inertia ().

Solution Set #1

Solution Set #2

subject to

J1    = 0.005;  % principal moment of inertia about x-axis
J2    = 0.012;  % principal moment of inertia about y-axis
J3    = 0.013;  % principal moment of inertia about z-axis
omega10 = 1;    % initial omega1
omega20 = 2;    % initial omega2
omega30 = 3;    % initial omega3
tau1  = - (J2 - J3)*omega20*omega30  % torque about x-axis
tau2  = - (J3 - J1)*omega30*omega10  % torque about y-axis
tau3  = - (J1 - J2)*omega10*omega20  % torque about z-axis

Verifying the solution (set #1):

omega1 = (sqrt((J2 - J3)*tau2*tau3))/(sqrt((J1 - J2)*(J1 - J3)*tau1))
omega2 = (tau3*sqrt((J1 - J2)*(J1 - J3)*tau1))/((J2 - J1)*sqrt((J2 - J3)*tau2*tau3))
omega3 = (tau2*sqrt((J1 - J2)*(J1 - J3)*tau1))/((J1 - J3)*sqrt((J2 - J3)*tau2*tau3))

omega1 =

1

omega2 =

2

omega3 =

3

Answer 2

[Paul]

Hi Vitor,

To make sure I understand the question, it assumes that a rigid body can have a constant, non-zero angular velocity vector resolved in body-fixed coordinates while being subjected to a non-zero torque, also resolved in body-fixed coordinates.

Let's try to solve this problem for the simple case where the body-fixed frame is defined by the principal axes

syms w [3 1] real
syms T [3 1] real
syms J1_1 J2_2 J3_3 positive
eq = T == cross(w,diag([J1_1 J2_2 J3_3])*w);
sol = solve(eq,w,'ReturnConditions',true,'Real',true)

sol = struct with fields:
            w1: [15×1 sym]
            w2: [15×1 sym]
            w3: [15×1 sym]
    parameters: [u    x    y    z]
    conditions: [15×1 sym]

So Matlab found 15 solutions under various conditions. Let's look at the first three conditions

sol.conditions(1:3)

ans = 

So the first three conditions all have the three torque components equal to zero, i..e., torque free motion. But those solutions are not of interest and so can be discarded.

Let's look at the last solution

sol.conditions(end)

ans = 

[sol.w1(end) sol.w2(end) sol.w3(end)]

ans = 

So that might be a solution of interest for some constant z for the case where T3 = 0 and J1_1 = J2_2 and J2_2 ~= J3_3.

Presumably you can go through all of the sol.conditions and find the ones of interest.

If you try to expand this approach to the general case, you'll need to enforce all of constraints on the inertia tensor, like symmetric and positive definite to begin with. There may be others that I'm forgetting.

Having said that, it seems to me that you don't need to solve the general case explicitly. Because the w and T vectors are constant in the body-fixed, principal axes frame, they are constant in any body-fixed frame. So start with whatever body-fixed frame you have, derive the direction cosine matrix from that frame to the principal axes, tranform the intertia tensor to principal axes, solve the problem in the principal axes frame, then tansform the solution back to the original body-fixed frame.

Interesting problem.