[MBDyn-users] How to compute the derivative of a rotation angle versus the rotation parameters (needed for the jacobian matrix in AssJac)

masarati at aero.polimi.it masarati at aero.polimi.it
Mon Feb 28 22:06:56 CET 2011


> Dear MBDyn users,
> Dear Mr. Masarati,
>
> I have the following problem:
> I need to impose a force at a specific node which depends on an rotation
> angle by means of an user defined element.
> Let's say Fg = f( Phi1 - Phi2 ).
> Where Fg is the force and Phi1 - Phi2 is the difference of the rotation
> angle between two nodes (e.g. Phi1 ... crank node, Phi2 ... stator node).
> The angles Phi1 and Phi2 are in relation to the reference frame of the
> stator node.
> In order to compute the difference between Phi1 and Phi2 I have written
> the following C++ code.
>
> // CrankAngle ... returns the angle ( delta_Phi = Phi1 - Phi2 ) between
> the 1-axis of the crank node and the 1-axis of the rotor node
> double unsteady_cylinder_gas_force::CrankAngle()const
> {
>             // matCrankRCurr ... orientation matrix of the crank in the
> reference frame of the stator
>             const Mat3x3 matCrankRCurr =
> m_pStatorNode->GetRCurr().Transpose() *
> m_pCrankNode->GetRCurr();
>             return atan2(matCrankRCurr(2,1),matCrankRCurr(1,1));
> }
>
> In this code it is assumed that the 3-axis of the crank node and the
> 3-axis of the stator node are nearly parallel but the spatial rotation of
> the stator node may be arbitrary large.
> Once the angle Phi1 - Phi2 is known, the force Fg can be evaluated and the
> residuum can be computed in AssRes.
> The problem is now how to compute the jacobian matrix in AssJac.
> Since the derivative Fgd of Fg versus the crank angle delta_Phi = Phi1 -
> Phi2 is known ( Fgd = diff(Fg(delta_Phi),delta_Phi) ) the inner derivative
> of delta_Phi versus the rotation parameters is needed in order to compute
> the jacobian matrix.
>
> The equation numbers of the rotation parameters are obtained by
> m_pCrankNode->iGetFirstPositionIndex()+4..6 and
> m_pStatorNode->iGetFirstPositionIndex()+4..6 where
> m_pCrankNode is a pointer to the structural node of the crank and
> m_pStatorNode is a pointer to the structural node of the stator.
>
>
> My approach was the following:
> I projected the rotation parameters to the reference frame of the stator
> because the 3-axis of the crank and the 3-axis of the rotor are nearly
> parallel as following:
> y^(stator) = diag(transpose(R_stator)) * y^(I)
> y^(stator) ... vector of the displacement and the rotation parameters of
> the crank node and the stator node in the reference frame of the stator
> y^(I)        ... vector of the displacement and the rotation parameters of
> the crank node and the stator node in the global reference frame
>      provided by the equation numbers
> m_pCrankNode->GetFristPositionIndex() and
> m_pStatorNode->GetFirstPositionIndex()
>      y^(I) = transpose( X_crank^(I)(1), X_crank^(I)(2), X_crank^(I)(3),
> g_crank^(I)(1), g_crank^(I)(2), g_crank^(I)(3), X_stator^(I)(1),
> X_stator^(I)(2), X_stator^(I)(3), g_stator^(I)(1), g_stator^(I)(2),
> g_stator^(I)(3) )
>      where X_crank^(I) and X_stator^(I) are the position of the crank node
> and stator node respectively in the global reference frame, and
> Phi_crank^(I) and Phi_stator^(I) are the rotation parameters of the
> crank node and the stator node in the global reference frame.
> R_stator ... rotation matrix of the stator node provided by
> m_pStatorNode->GetRCurr()
> diag(Q) ... a function that returns a matrix with submatrix A in its main
> diagonal
>                                          | Q          |
>                 e.g. diag(Q) =  |     Q     |
>                                       |         Q |
>
> If the rotation parameters Phi_crank(3) and Phi_stator(3) were just
> angles, the angle delta_Phi = Phi1 - Phi2 can be written as
> delta_Phi = y^(stator)(6) -y^(stator)(12) = g_crank^(stator)(3) -
> g_stator^(stator)(3)
> and
> delta_Phi = ( 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, -1 ) * y^(stator) = A *
> diag(transpose(R_stator)) * y^(I)
> with A = ( 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, -1 ).
>
> Because you wrote in your mail from 25.01.2011 that the rotation
> parameters g_crank^(I) and g_stator^(I) are not just angles, my question
> is now how to compute the angle delta_Phi from the rotation parameters
> g_crank^(I) and g_stator^(I).
>
> For example if the rotation g parameters were Gibbs-Rodriguez parameters
> of the form
> g = 2*tan(Phi/2) as described in "Comprehensive Multibody AeroServoElastic
> Analysis of Integrated Rotorcraft Active Controls" chapter 3.2.2
> the angle delta_Phi could be computed as
> delta_Phi = 2 * ( atan(g_crank^(stator)(3)/2) -
> atan(g_stator^(stator)(3)/2 ).
>
> But in the C++ code of MBDyn it seems that the updated rotation matrices
> of structural nodes are computed by means of the method
> Mat3x3_Manip::Manipulate which is virtual and has different
> implementations.
> So it is not obvious for me how to compute the inner derivative of the
> angle delta_Phi versus the rotation parameters g.

Hello.  Please forgive me if I don't go through all the details of your
analysis.  Please find a simple (and brief) expression of the virtual
perturbation of the relative rotation of two nodes about a fixed axis. 
This is consistent with the way rotations are handled in MBDyn.  Please
let us know if there's anything unclear.

Cheers, p.
-------------- next part --------------
A non-text attachment was scrubbed...
Name: x.pdf
Type: application/pdf
Size: 19350 bytes
Desc: not available
URL: <http://mail.mbdyn.org/pipermail/mbdyn-users/attachments/20110228/72e3a49f/attachment-0001.pdf>


More information about the MBDyn-users mailing list