[MBDyn-users] How to compute the derivative of a rotation angle versus the rotation parameters (needed for the jacobian matrix in AssJac)
Resch Reinhard
Reinhard.RESCH at accomp.it
Mon Feb 28 11:00:24 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.
Cheers
Reinhard Resch
=========================
Research & Development
ACC Austria GesmbH
Jahnstrasse 30
8280 Fürstenfeld, Austria
Tel: ++43 3382 5010-434
Fax: ++43 3382 5010-355
www: http://www.the-acc-group.it/
mail: reinhard.resch at accomp.it
________________________________
Le informazioni contenute in questo messaggio e negli eventuali allegati sono riservate e per uso esclusivo del destinatario. Persone diverse dallo stesso non possono copiare o distribuire il messaggio a terzi. Chiunque riceva questo messaggio per errore, è pregato di distruggerlo e di informare immediatamente il mittente.
HINWEIS: Dies ist eine vertrauliche Nachricht und nur für den Adressaten bestimmt. Es ist nicht erlaubt, diese Nachricht zu kopieren oder Dritten zugänglich zu machen. Sollten Sie diese Nachricht irrtümlich erhalten haben, bitte ich um Ihre Mitteilung per E-Mail oder unter der oben angegebenen Telefonnummer.
The information transmitted is intended for the person or entity to which it is addressed and may contain confidential and/or privileged material. Any review, retransmission, dissemination or other use of, or taking of any action in reliance upon, this information by persons or entities other than the intended recipient is prohibited. If you received this in error, please contact the sender and delete the material from any computer.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.mbdyn.org/pipermail/mbdyn-users/attachments/20110228/9c8eb11e/attachment.htm>
More information about the MBDyn-users
mailing list