I) A Lagrangian variational principle for Euler's equations for a rigid body
$$ \tag{1} (DL)_i ~=~M_i, \qquad i\in\{1,2,3\}, $$
is e.g. explained in Ref. 1. Here the angular momentum $L_i$, $i\in\{1,2,3\}$, along the three principal axes of inertia is tied to the angular velocity $\omega_i$, $i\in\{1,2,3\}$, by the formula
$$\tag{2} L_i~:=~I_i \omega_i, \qquad i\in\{1,2,3\}, \qquad (\text{no sum over }i).$$
The covariant time-derivative $D$ of a vector $\eta_i$, $i\in\{1,2,3\}$, is defined as
$$\tag{3} (D\eta)_i ~:=~ \dot{\eta}_i+(\omega\times\eta)_i, \qquad i\in\{1,2,3\}. $$
The angular velocity vector $\omega$ plays the role of a non-Abelian gauge connection/potential.
II) To see the $so(3)$ Lie algebra, we map an infinitesimal rotation vector $\alpha$ into an antisymmetric real $3\times3$ matrix $r(\alpha)\in so(3)$ as
$$\tag{4} \alpha_i \quad\longrightarrow\quad r(\alpha)_{jk}~:=~\sum_{i=1}^3\alpha_i \varepsilon_{ijk}. $$
The $so(3)$ Lie-bracket is given by (minus) the vector cross product
$$\tag{5} [r(\alpha),r(\beta)]~=~r(\beta\times \alpha). $$
Similarly, for the corresponding $SO(3)$ Lie group, a finite rotation vector $\alpha$ maps into an orthogonal $3\times3$ rotation matrix $R(\alpha)\in SO(3)$ as explained in this Phys.SE post. Infinitesimally, for an infinitesimal rotation $|\delta\alpha| \ll 1$, the correspondence is
$$\tag{6} R(\delta\alpha)_{jk} ~=~\delta_{jk} + r(\delta\alpha)_{jk} + {\cal O}(\delta\alpha^2). $$
III) A finite non-Abelian gauge transformation
$\omega\longrightarrow\omega^{\alpha}$ takes the form
$$\tag{7} r(\omega^{\alpha})~=~R(-\alpha) \left(\frac{d}{dt}-r(\omega)\right)R(\alpha), \qquad \alpha\in \mathbb{R}^3.$$
An infinitesimal non-Abelian gauge transformation $\delta$ takes the form
$$\tag{8} r(\delta\omega)~=~\frac{d}{dt}r(\delta\alpha)-[r(\omega),r(\delta\alpha)],$$
or equivalently
$$\tag{9} \delta\omega_i~=~(D\delta\alpha)_i, \qquad i\in\{1,2,3\}, $$
where $\delta\alpha$ denotes an infinitesimal rotation vector corresponding to an $so(3)$ Lie algebra element $r(\delta\alpha)$.
We call (7)-(9) gauge transformations for semantic reasons, because of their familiar form, but note that (most of) them are not unphysical/spurious transformations. We stress that the angular velocity $\omega$ is a physical variable.
IV) Finally we are ready to discuss the action principle. The finite rotation vector $\alpha(t)\in \mathbb{R}^3$ plays the role of independent dynamical variables for the action principle. One may think of virtual rotation paths $\alpha:[t_i,t_f]\to \mathbb{R}^3$ as a reparametrization of virtual angular velocity paths $\omega:[t_i,t_f]\to \mathbb{R}^3$. The action reads
$$\tag{10} S[\alpha,\omega]~=~\int_{t_i}^{t_f} \! dt ~L $$
with Lagrangian
$$\tag{11} L~=~\frac{1}{2} L^{\alpha}\cdot \omega^{\alpha} + M\cdot \alpha , $$
where
$$\tag{12} L^{\alpha}_i~:=~I_i \omega^{\alpha}_i, \qquad i\in\{1,2,3\}, \qquad (\text{no sum over }i).$$
The Lagrangian (11) consists of rotational kinetic energy plus a source term from the torque $M$. Here $\omega^{\alpha}$ is the actual angular velocity vector, while $\omega$ here (in contrast to above) is a fixed non-dynamical reference vector, which is not varied. It is sort of a gauge-fixing choice. Infinitesimal variation yields
$$ \tag{13}\delta L
~\stackrel{(11)}{=}~ L^{\alpha}\cdot \delta\omega^{\alpha}
+ M\cdot \delta\alpha
~\stackrel{(9)}{=}~ L^{\alpha}\cdot \left(\frac{d}{dt}\delta\alpha
+ (\omega^{\alpha}\times\delta\alpha)\right)
+ M\cdot \delta\alpha ,$$
which (after integration by parts and appropriate boundary conditions) leads to Euler's equations (1) for the angular velocity vector $\omega^{\alpha}$.
References:
- J.E. Marsden and T.S. Ratiu, Introduction to
Mechanics and Symmetry, 1998.
I did not follow your derivation of $\frac{{\rm d}\mathbf{I}}{{\rm d}t}$. In most textbooks it is evaluated as follows $$\frac{{\rm d}\mathbf{I}}{{\rm d}t} =\boldsymbol{ \omega } \times \mathbf{I} = \begin{vmatrix} 0 & -\omega_z & \omega_y \\ \omega_z & 0 & -\omega_x \\ -\omega_y & \omega_x & 0 \end{vmatrix} \begin{vmatrix} I_{xx} & I_{xy} & I_{xz} \\ I_{xy} & I_{yy} & I_{yz} \\ I_{xz} & I_{yz} & I_{zz} \end{vmatrix} $$
with the added caveat that $\mathbf{I}$ depends on the orientation of the body. The orientation may be tracked using Euler angles, Quaternions or just the 3×3 rotation matrix $\mathbf{R}$. Either way the end result is that the mass moment of inertia tensor needs to be calculated at every instant from the MMOI in the body frame
$$ \mathbf{I} = \mathbf{R}\,\mathbf{I}_{\rm body} \,\mathbf{R}^\top $$
In the end you have the equations of motion
$$ \left. \boldsymbol{\tau} = \mathbf{I}\, \boldsymbol{\dot{\omega}} + \boldsymbol{\omega} \times \mathbf{I} \boldsymbol{\omega}\;\; \right\} \;\; \boldsymbol{\dot{\omega}} = \mathbf{I}^{-1}\left(\boldsymbol{\tau} - \boldsymbol{\omega} \times \mathbf{I} \boldsymbol{\omega} \right) $$
It is also common to express the above in terms of angular momentum in the following algorithm. Each integration step is given the rotation matrix $\mathbf{R}$ and momentum vector $\boldsymbol{L}$
$$
\begin{array}{c|cc}
\text{Step} & \text{Calculation} & \text{Notes}\\
\hline 0 & \mathbf{I}=\mathbf{R}\mathbf{I}_{{\rm body}}\mathbf{R}^{\top} & \text{MMOI in world coorinates}\\
1 & \boldsymbol{\omega}=\mathbf{I}^{-1}\boldsymbol{L} & \text{Extract rotational vector}\\
2 & \dot{\mathbf{R}}=\boldsymbol{\omega}\times\mathbf{R} & \text{Change in rotation}^\star\\
3 & \dot{\boldsymbol{L}}=\boldsymbol{\tau}(t,\mathbf{R},\boldsymbol{\omega}) & \text{Change in momentum due to torque }\boldsymbol{\tau}
\end{array}
$$
* Note: When integrating the rotation matrix $\mathbf{R}$ using Runge-Kutta the result of $\mathbf{R} \rightarrow \mathbf{R} + h \dot{\mathbf{R}}$ is no longer a rotation matrix and the solution will quickly degrade in accuracy.
So instead, people often use quaternions $\boldsymbol{\hat{q}} = \pmatrix{ \boldsymbol{q}_{\rm v} & q_{\rm s}} $ that describe the rotation as $$ \mathbf{R} = \mathbf{1} + 2 q_{\rm s} [ \boldsymbol{q}_{\rm v}\times] + 2 [ \boldsymbol{q}_{\rm v} \times][ \boldsymbol{q}_{\rm v} \times] $$ where $[ \boldsymbol{q}_{\rm v} \times] = \begin{vmatrix} 0 & -z & y \\ z & 0 & -x \\ -y & x & 0 \end{vmatrix}$ is the 3×3 cross product matrix operator of the vector part of the quaternion $\boldsymbol{q}_{\rm v}$.
The derivative of the quaternion is defined as $$ \dot{\boldsymbol{\hat{q}}} = \frac{1}{2} \pmatrix{ -\boldsymbol{\omega}^\top \boldsymbol{q}_{\rm v} \\ q_{\rm s} \boldsymbol{\omega} + \boldsymbol{\omega} \times \boldsymbol{q}_{\rm v} }$$
But often still people get this step wrong, because integrating the above $\boldsymbol{\hat{q}} \rightarrow \boldsymbol{\hat{q}} + h \dot{\boldsymbol{\hat{q}}}$ still brakes the rotation representation.
The proper way to take an integration step with quaternions is as follows. Given $\boldsymbol{\hat{q}} = \pmatrix{\boldsymbol{q}_{\rm v} & q_{\rm s}}$ and $\boldsymbol{\omega}$ vector
$$ \pmatrix{\boldsymbol{q}_{\rm v} \\ q_{\rm s}} \rightarrow \begin{vmatrix}
\cos(\tfrac{\theta}{2} ) & -\sin(\tfrac{\theta}{2} ) \boldsymbol{z}^\top \\
\sin(\tfrac{\theta}{2} ) \boldsymbol{z} & \cos(\tfrac{\theta}{2} ) + \sin(\tfrac{\theta}{2} ) [\boldsymbol{z}\times]
\end{vmatrix} \pmatrix{\boldsymbol{q}_{\rm v} \\ q_{\rm s}} $$
where $\theta = h \| \boldsymbol{\omega} \|$ is the step angle and $\boldsymbol{z} = \boldsymbol{\omega}/\|\boldsymbol{\omega}\|$ is the step rotation axis.
The resulting quaternion still represents rotations always and does not drift away like other formulations do.
Best Answer
tl;dr: Yes, this is valid. Getting there is nontrivial.
The $\boldsymbol \omega_B \times (\mathrm I_B \, \boldsymbol \omega_B)$ term in $\boldsymbol \tau_B = \mathrm I_B \, \dot{\boldsymbol \omega}_B + \boldsymbol \omega_B \times (\mathrm I_B \, \boldsymbol \omega_B)$ results from working in a rotating frame. It is a fictitious torque, very much akin to the fictitious forces that arise with Newton's second law when computed in a rotating frame. Just as those fictitious forces vanish in an inertial frame, so does the fictitious torque. The inertial frame rotational equation of motion is
$$\boldsymbol\tau_I = \frac d{dt} (\mathrm I_I\, \boldsymbol \omega_I) = \mathrm I_I \, \dot {\boldsymbol \omega}_I + \dot {\mathrm I}_I \,\boldsymbol \omega_I$$
This raises the question, what is the time derivative of the inertia tensor from the perspective of an inertial frame? To determine this, the relationship between the inertia tensor as expressed in the body frame and as expressed inertial frame is needed. This is $$\mathrm I_I = \mathrm T_{R\to I}\,\mathrm I_B\,\mathrm T_{R\to I}^{\,T}$$ where $\mathrm I_I$ is the inertia tensor expressed in the inertial frame, $\mathrm T_{R\to I}$ is the transformation matrix from the body frame to the inertial frame, $\mathrm I_B$ the inertia tensor expressed in the body frame, and $\mathrm T_{R\to I}^{\,T}$ is the transpose of (and also the inverse of) $\mathrm T_{R\to I}$. Taking the derivative with respect to time yields $$\dot {\mathrm I}_I = \dot {\mathrm T}_{R\to I} \, \mathrm I_B\,\mathrm T_{R\to I}^{\,T} + \mathrm T_{R\to I} \, \mathrm I_B \, \dot {\mathrm T}_{R\to I}^{\,T}$$ Without derivation, the time derivative of the transformation matrix is $$\dot {\mathrm T}_{R\to I} = \mathrm T_{R\to I} \operatorname{Sk}(\boldsymbol \omega_B) = \operatorname{Sk}(\boldsymbol \omega_I)\mathrm T_{R\to I}$$ where $\operatorname{Sk}(a)$ is the skew symmetric cross product matrix generated by the vector $\boldsymbol a$: $$\operatorname{Sk}(\boldsymbol a) \equiv \begin{pmatrix} \phantom{-}0 & -a_z & \phantom{-}a_y \\ \phantom{-}a_z & \phantom{-}0 & -a_x \\ -a_y & \phantom{-}a_x & \phantom{-}0 \end{pmatrix}$$ Note that $\operatorname{Sk}(\boldsymbol a) \, \boldsymbol b = \boldsymbol a \times \boldsymbol b$. With the above, the time derivative of the inertia tensor as expressed in the inertial frame is $$\dot {\mathrm I}_I = \operatorname{Sk}(\boldsymbol \omega_I) \mathrm T_{R\to I} \, \mathrm I_B\,\mathrm T_{R\to I}^{\,T} + \mathrm T_{R\to I} \, \mathrm I_B \, \mathrm T_{R\to I}^{\,T} \operatorname{Sk}^T(\boldsymbol \omega_I) = \operatorname{Sk}(\boldsymbol \omega_I)\,\mathrm I_I - \mathrm I_I \, \operatorname{Sk}(\boldsymbol \omega_I)$$ This can be nonzero due to the non-commutative nature of matrix multiplication. We're not interested so much in the form of $\dot {\mathrm I}_I$ in general so much as the product $\dot {\mathrm I}_I\,\boldsymbol \omega_I$. This is $$\dot {\mathrm I}_I\,\boldsymbol \omega_I = \operatorname{Sk}(\boldsymbol \omega_I)\,\mathrm I_I \,\boldsymbol \omega_I - \mathrm I_I \, \operatorname{Sk}(\boldsymbol \omega_I) \,\boldsymbol \omega_I = \boldsymbol \omega_I \times (\mathrm I_I \,\boldsymbol \omega_I) - \mathrm I_I \, (\boldsymbol \omega_I \times \boldsymbol \omega_I) = \boldsymbol \omega_I \times (\mathrm I_I \,\boldsymbol \omega_I)$$ and thus $$\boldsymbol\tau_I = \mathrm I_I \, \dot {\boldsymbol \omega}_I + \boldsymbol \omega_I \times (\mathrm I_I \,\boldsymbol \omega_I)$$