[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() *
> 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
> 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
> diag(Q) ... a function that returns a matrix with submatrix A in its main
> | 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) -
> 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
> 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.
-------------- next part --------------
A non-text attachment was scrubbed...
Size: 19350 bytes
Desc: not available
More information about the MBDyn-users