MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
module-asynchronous_machine.cc
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/modules/module-asynchronous_machine/module-asynchronous_machine.cc,v 1.11 2017/01/12 14:47:25 masarati Exp $ */
2 /*
3  * MBDyn (C) is a multibody analysis code.
4  * http://www.mbdyn.org
5  *
6  * Copyright (C) 1996-2017
7  *
8  * Pierangelo Masarati <masarati@aero.polimi.it>
9  * Paolo Mantegazza <mantegazza@aero.polimi.it>
10  *
11  * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
12  * via La Masa, 34 - 20156 Milano, Italy
13  * http://www.aero.polimi.it
14  *
15  * Changing this copyright notice is forbidden.
16  *
17  * This program is free software; you can redistribute it and/or modify
18  * it under the terms of the GNU General Public License as published by
19  * the Free Software Foundation (version 2 of the License).
20  *
21  *
22  * This program is distributed in the hope that it will be useful,
23  * but WITHOUT ANY WARRANTY; without even the implied warranty of
24  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25  * GNU General Public License for more details.
26  *
27  * You should have received a copy of the GNU General Public License
28  * along with this program; if not, write to the Free Software
29  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
30  */
31 /*
32  AUTHOR: Reinhard Resch <r.resch@secop.com>
33  Copyright (C) 2011(-2017) all rights reserved.
34 
35  The copyright of this code is transferred
36  to Pierangelo Masarati and Paolo Mantegazza
37  for use in the software MBDyn as described
38  in the GNU Public License version 2.1
39 */
40 
70 #include "mbconfig.h" // This goes first in every *.c,*.cc file
71 
72 #include <cassert>
73 #include <cstdio>
74 #include <cmath>
75 #include <cfloat>
76 #include <iostream>
77 #include <iomanip>
78 #include <limits>
79 
80 #include "dataman.h"
81 #include "userelem.h"
83 
85 : virtual public Elem, public UserDefinedElem
86 {
87 public:
88  asynchronous_machine(unsigned uLabel, const DofOwner *pDO,
89  DataManager* pDM, MBDynParser& HP);
90  virtual ~asynchronous_machine(void);
91  virtual void Output(OutputHandler& OH) const;
92  virtual unsigned int iGetNumDof(void) const;
93  virtual DofOrder::Order GetDofType(unsigned int i) const;
94  virtual std::ostream& DescribeDof(std::ostream& out, const char *prefix, bool bInitial) const;
95  virtual std::ostream& DescribeEq(std::ostream& out, const char *prefix, bool bInitial) const;
96  virtual unsigned int iGetNumPrivData(void) const;
97  virtual unsigned int iGetPrivDataIdx(const char *s) const;
98  virtual doublereal dGetPrivData(unsigned int i) const;
99  virtual void SetInitialValue(VectorHandler& X);
100  virtual void Update(const VectorHandler& XCurr,const VectorHandler& XPrimeCurr);
101  virtual void AfterPredict(VectorHandler& X, VectorHandler& XP);
102  virtual void WorkSpaceDim(integer* piNumRows, integer* piNumCols) const;
105  doublereal dCoef,
106  const VectorHandler& XCurr,
107  const VectorHandler& XPrimeCurr);
109  AssRes(SubVectorHandler& WorkVec,
110  doublereal dCoef,
111  const VectorHandler& XCurr,
112  const VectorHandler& XPrimeCurr);
113  virtual int iGetNumConnectedNodes(void) const;
114  virtual void GetConnectedNodes(std::vector<const Node *>& connectedNodes) const;
115  virtual void SetValue(DataManager *pDM, VectorHandler& X, VectorHandler& XP,
117  virtual std::ostream& Restart(std::ostream& out) const;
118  virtual unsigned int iGetInitialNumDof(void) const;
119  virtual void
120  InitialWorkSpaceDim(integer* piNumRows, integer* piNumCols) const;
121  virtual VariableSubMatrixHandler&
123  const VectorHandler& XCurr);
124  virtual SubVectorHandler&
125  InitialAssRes(SubVectorHandler& WorkVec, const VectorHandler& XCurr);
126 
127 private:
128  bool IsMotorOn(void) const;
129 
130 private:
142  static const doublereal sm_SingTol;
143 };
144 
145 const doublereal asynchronous_machine::sm_SingTol = std::pow(std::numeric_limits<doublereal>::epsilon(), 0.9);
146 
158  unsigned uLabel, const DofOwner *pDO,
159  DataManager* pDM, MBDynParser& HP)
160 : Elem(uLabel, flag(0)),
161  UserDefinedElem(uLabel, pDO),
162  m_pRotorNode(0),
163  m_pStatorNode(0),
164  m_MK(0.),
165  m_sK(0.),
166  m_OmegaS(0),
167  m_MotorOn(0),
168  m_M(0.),
169  m_dM_dt(0.),
170  m_dM_dt2(0.),
171  m_omega(0.),
172  m_domega_dt(0.)
173 {
174  if (HP.IsKeyWord("help")) {
175  silent_cout(
176  "\n"
177  "Module: asynchronous_machine\n"
178  "\n"
179  " This module implements an asynchronous electric machine\n"
180  " according to\n"
181  "\n"
182  " Maschinendynamik\n"
183  " Hans Dresig, Franz Holzwei├čig\n"
184  " Springer, 2009\n"
185  " page 58 equation 1.122 and 1.123\n"
186  "\n"
187  "Syntax:\n"
188  " asynchronous_machine,\n"
189  " rotor, (label) <rotor_node>,\n"
190  " stator, (label) <stator_node>,\n"
191  " MK, (real) <MK>,\n"
192  " sK, (real) <sK>,\n"
193  " OmegaS, (DriveCaller) <omega_s>\n"
194  " [ , motor on , (DriveCaller) <motor_on> ]\n"
195  " [ , M0 , (real) <M0> ]\n"
196  " [ , MP0 , (real) <MP0> ]\n"
197  "\n"
198  " MK ... breakdown torque ( MK > 0 )\n"
199  " sK ... breakdown slip ( sK > 0 )\n"
200  " OmegaS = 2 * pi * f / p * sense_of_rotation ... synchronous angular velocity\n"
201  " motor_on ... power supply is turned off when 0, on otherwise\n"
202  " M0 ... initial torque ( default 0 )\n"
203  " MP0 ... initial torque derivative ( default 0 )\n"
204  "\n"
205  " The axis of rotation is assumed to be axis 3 of the reference frame of the stator node.\n"
206  << std::endl);
207 
208  if (!HP.IsArg()) {
209  // Exit quietly if nothing else is provided
210  throw NoErr(MBDYN_EXCEPT_ARGS);
211  }
212  }
213 
214  if (!HP.IsKeyWord("rotor")) {
215  silent_cerr("asynchronous_machine(" << GetLabel() << "): keyword \"rotor\" expected at line " << HP.GetLineData() << std::endl);
217  }
218 
220 
221  if (!m_pRotorNode) {
222  silent_cerr("asynchronous_machine(" << GetLabel() << "): structural node expected at line " << HP.GetLineData() << std::endl);
224  }
225 
226  if (!HP.IsKeyWord("stator")) {
227  silent_cerr("asynchronous_machine(" << GetLabel() << "): keyword \"stator\" expected at line " << HP.GetLineData() << std::endl);
229  }
230 
232 
233  if (!m_pStatorNode) {
234  silent_cerr("asynchronous_machine(" << GetLabel() << "): structural node expected at line " << HP.GetLineData() << std::endl);
236  }
237 
238  if (!HP.IsKeyWord("MK")) {
239  silent_cerr("asynchronous_machine(" << GetLabel() << "): keyword \"MK\" expected at line " << HP.GetLineData() << std::endl);
241  }
242  m_MK = HP.GetReal();
243 
244  if (!HP.IsKeyWord("sK")) {
245  silent_cerr("asynchronous_machine(" << GetLabel() << "): keyword \"sK\" expected at line " << HP.GetLineData() << std::endl);
247  }
248  m_sK = HP.GetReal();
249 
250  if (!HP.IsKeyWord("OmegaS")) {
251  silent_cerr("asynchronous_machine(" << GetLabel() << "): keyword \"OmegaS\" expected" << std::endl);
253  }
255 
256  if (HP.IsKeyWord("motor" "on") || HP.IsKeyWord("motor_on")) {
258  } else {
260  }
261 
262  if (HP.IsKeyWord("M0")) {
263  m_M = HP.GetReal();
264  }
265 
266  if (HP.IsKeyWord("MP0")) {
267  m_dM_dt = HP.GetReal();
268  }
269 
271 
272  Vec3 e3(m_pStatorNode->GetRCurr().GetCol(3));
273  const Vec3& omega1 = m_pRotorNode->GetWCurr();
274  const Vec3& omega2 = m_pStatorNode->GetWCurr();
275 
276  m_omega = e3.Dot(omega1 - omega2);
277 
278  const doublereal OmegaS = m_OmegaS.dGet();
279 
280  pDM->GetLogFile() << "asynchronous machine: "
281  << uLabel << " "
282  << m_pRotorNode->GetLabel() << " "
283  << m_pStatorNode->GetLabel() << " "
284  << m_MK << " "
285  << m_sK << " "
286  << OmegaS << " "
287  << IsMotorOn() << " "
288  << m_M << " "
289  << m_dM_dt << " "
290  << m_omega
291  << std::endl;
292 
293  // Note: The ODE for the asynchronous machine is defined for positive sign of OmegaS only!
294  // At user level the signs of M, dM/dt and d^2M/dt^2 are defined to be positive
295  // in positive coordinate system direction in the reference frame of the stator.
296  m_M *= copysign(1., OmegaS);
297  m_dM_dt *= copysign(1., OmegaS);
298  m_dM_dt2 *= copysign(1., OmegaS);
299 }
300 
310 void
312 {
313  const integer iFirstIndex = iGetFirstIndex();
314 
315  const integer intTorqueDerivativeRowIndex = iFirstIndex + 1;
316  const integer intTorqueRowIndex = iFirstIndex + 2;
317  const integer intOmegaRowIndex = iFirstIndex + 3;
318 
319  // save the current state needed by dGetPrivData() and Output()
320  m_M = X(intTorqueRowIndex);
321  m_dM_dt = X(intTorqueDerivativeRowIndex);
322  m_dM_dt2 = XP(intTorqueDerivativeRowIndex);
323  m_omega = X(intOmegaRowIndex);
324  m_domega_dt = XP(intOmegaRowIndex);
325 }
326 
328 {
329  Update(X, XP);
330 }
331 
337 bool
339 {
340  return (m_MotorOn.dGet() != 0.);
341 }
342 
344 {
345  // destroy private data
346 #ifdef DEBUG
347  std::cerr << __FILE__ << ":" << __LINE__ << ":" << __PRETTY_FUNCTION__ << std::endl;
348 #endif
349 }
350 
369 void
371 {
372  if (bToBeOutput()) {
373  // Note: The ODE for the asynchronous machine is defined for positive sign of OmegaS only!
374  // At user level the signs of M, dM/dt and d^2M/dt^2 are defined to be positive
375  // in positive coordinate system direction in the reference frame of the stator.
376  doublereal OmegaS = m_OmegaS.dGet();
377  doublereal M = m_M*copysign(1., OmegaS);
378  doublereal dM_dt = m_dM_dt*copysign(1., OmegaS);
379  doublereal dM_dt2 = m_dM_dt2*copysign(1., OmegaS);
380 
382  // output the current state: the column layout is as follows
383 
384  // 1 2 3 4 5 6 7
385  // M dM_dt dM_dt2 omega domega_dt OmegaS gamma
386  //
387  // 0 label of the element
388  // 1 motor torque
389  // 2 derivative of the motor torque versus time
390  // 3 second derivative of the motor torque versus time
391  // 4 angular velocity of the rotor node relative to the stator node around axis 3
392  // 5 angular acceleration of the rotor node relative to the stator node around axis 3
393  // 6 synchronous angular velocity (might be a function of the time in case of an frequency inverter)
394  // 7 power supply is turned on or off
395 
396  OH.Loadable() << GetLabel()
397  << " " << M
398  << " " << dM_dt
399  << " " << dM_dt2
400  << " " << m_omega
401  << " " << m_domega_dt
402  << " " << OmegaS
403  << " " << IsMotorOn()
404  << std::endl;
405  }
406  }
407 }
408 
414 void
416 {
417  *piNumRows = 9;
418  *piNumCols = 9;
419 }
420 
425 unsigned int
427 {
428  return 3;
429 }
430 
437 asynchronous_machine::GetDofType(unsigned int i) const
438 {
439 
440  return DofOrder::DIFFERENTIAL;
441 }
442 
451 std::ostream&
452 asynchronous_machine::DescribeDof(std::ostream& out, const char *prefix, bool bInitial) const
453 {
454  integer iIndex = iGetFirstIndex();
455 
456  out
457  << prefix << iIndex + 1 << ": motor torque derivative [MP]" << std::endl
458  << prefix << iIndex + 2 << ": motor torque [M]" << std::endl
459  << prefix << iIndex + 3 << ": relative angular velocity [omega]" << std::endl;
460 
461  return out;
462 }
463 
472 std::ostream&
473 asynchronous_machine::DescribeEq(std::ostream& out, const char *prefix, bool bInitial) const
474 {
475  integer iIndex = iGetFirstIndex();
476 
477  out
478  << prefix << iIndex + 1 << ": motor DAE [f1]" << std::endl
479  << prefix << iIndex + 2 << ": motor torque derivative [f2]" << std::endl
480  << prefix << iIndex + 3 << ": angular velocity derivative [f9]" << std::endl;
481 
482  return out;
483 }
484 
491 unsigned int
493 {
494  return 5;
495 }
496 
511 unsigned int
513 {
514  static const struct {
515  int index;
516  char name[7];
517  }
518 
519  data[] = {
520  { 1, "M" }, // motor torque
521  { 2, "MP"}, // derivative of the motor torque versus time diff(M,t) (named MP instead of dM_dt for compatibility with other elements)
522  { 3, "MPP"}, // second derivative of the motor torque versus time diff(M,t,2) (named MPP instead of dM_dt2 for compatibility with other elements)
523  { 4, "omega"}, // angular velocity of the rotor node relative to the stator node
524  { 5, "omegaP"} // angular acceleration of the rotor node relative to the stator node (named omegaP instead of domega_dt for compatibility with other elements)
525  };
526 
527  for (unsigned i = 0; i < sizeof(data) / sizeof(data[0]); ++i ) {
528  if (0 == strcmp(data[i].name,s)) {
529  return data[i].index;
530  }
531  }
532 
533  silent_cerr("asynchronous_machine(" << GetLabel() << "): no private data \"" << s << "\"" << std::endl);
534 
535  return 0;
536 }
537 
545 {
546  // Note: The ODE for the asynchronous machine is defined for positive sign of OmegaS only!
547  // At user level the signs of M, dM/dt and d^2M/dt^2 are defined to be positive
548  // in positive coordinate system direction in the reference frame of the stator.
549 
550  const doublereal OmegaS = m_OmegaS.dGet();
551 
552  switch (i) {
553  case 1:
554  return m_M*copysign(1., OmegaS);
555  case 2:
556  return m_dM_dt*copysign(1., OmegaS);
557  case 3:
558  return m_dM_dt2*copysign(1., OmegaS);
559  case 4:
560  return m_omega;
561  case 5:
562  return m_domega_dt;
563  }
564 
566 }
567 
571 void
573 {
574  return;
575 }
576 
628  doublereal dCoef,
629  const VectorHandler& XCurr,
630  const VectorHandler& XPrimeCurr)
631 {
632 #ifdef DEBUG
633  std::cerr << __FILE__ << ':' << __LINE__ << ':' << __PRETTY_FUNCTION__ << std::endl;
634 #endif
635 
636  integer iNumRows = 0;
637  integer iNumCols = 0;
638 
639  WorkSpaceDim(&iNumRows, &iNumCols);
640 
641  FullSubMatrixHandler& WorkMat = WorkMat_.SetFull();
642 
643  WorkMat.ResizeReset(iNumRows, iNumCols);
644 
645  const integer iFirstIndex = iGetFirstIndex();
646 
647  const integer intTorqueDerivativeColumnIndex = iFirstIndex + 1;
648  const integer intTorqueColumnIndex = iFirstIndex + 2;
649  const integer intOmegaColumnIndex = iFirstIndex + 3;
650 
651  const integer intTorqueDerivativeRowIndex = iFirstIndex + 1;
652  const integer intTorqueRowIndex = iFirstIndex + 2;
653  const integer intOmegaRowIndex = iFirstIndex + 3;
654 
655  const integer intRotorPositionIndex = m_pRotorNode->iGetFirstPositionIndex();
656  const integer intStatorPositionIndex = m_pStatorNode->iGetFirstPositionIndex();
657 
658  const integer intRotorMomentumIndex = m_pRotorNode->iGetFirstMomentumIndex();
659  const integer intStatorMomentumIndex = m_pStatorNode->iGetFirstMomentumIndex();
660 
661  /*
662  * M ... motor torque
663  * diff(M,t) ... derivative of the motor torque versus time
664  * g1 ... rotation parameter of the rotor
665  * g2 ... rotation parameter of the stator
666  * R2 ... rotation matrix of the stator node
667  * omega1 ... rotor angular velocity
668  * omega2 ... stator angular velocity
669  *
670  * e3 = ( 0, 0, 1 )^T ... axis of rotation in the reference frame of the stator node
671  *
672  * omega = e3^T * R2^T * ( omega1 - omega2 )
673  *
674  * s = 1 - omega / OmegaS
675  *
676  * | y1 | | diff(M,t) |
677  * | y2 | | M |
678  * y = | y3 | = | g1 |
679  * | y4 | | g2 |
680  * | y5 | | omega |
681  *
682  * | f1 | | diff(y1,t) + ( 2 * sK + sign(OmegaS) * diff(y5,t) / ( s * OmegaS^2 ) ) * y1 * abs(OmegaS) + ( ( sK^2 + s^2 ) * OmegaS^2 + sign(OmegaS) * diff(y5,t) * sK / s ) * y2 - 2 * MK * sK * s * OmegaS^2 |
683  * | f2 | | diff(y2,t) - y1 | |
684  * f = | f3 | = | R2 * e3 * y2 |
685  * | f4 | | -R2 * e3 * y2 |
686  * | f5 | | y5 - e3^T * R2^T * ( omega1 - omega2 ) |
687  *
688  * 1, 2, 3, 6, 9
689  * | df1/dy1, df1/dy2, df1/dy3, df1/dy4, df1/dy5 | 1
690  * | df2/dy1, df2/dy2, df2/dy3, df2/dy4, df2/dy5 | 2
691  * diff(f,y) = | df3/dy1, df3/dy2, df3/dy3, df3/dy4, df3/dy5 | 3
692  * | df4/dy1, df4/dy2, df4/dy3, df4/dy4, df4/dy5 | 6
693  * | df5/dy1, df5/dy2, df5/dy3, df5/dy4, df5/dy5 | 9
694  *
695  * WorkMat = -diff(f,diff(y,t)) - dCoef * diff(f,y)
696  */
697  WorkMat.PutRowIndex(1, intTorqueDerivativeRowIndex);
698  WorkMat.PutColIndex(1, intTorqueDerivativeColumnIndex);
699 
700  WorkMat.PutRowIndex(2, intTorqueRowIndex);
701  WorkMat.PutColIndex(2, intTorqueColumnIndex);
702 
703  for (int iCnt = 1; iCnt <= 3; ++iCnt) {
704  WorkMat.PutRowIndex(iCnt + 2, intRotorMomentumIndex + iCnt + 3);
705  WorkMat.PutColIndex(iCnt + 2, intRotorPositionIndex + iCnt + 3);
706 
707  WorkMat.PutRowIndex(iCnt + 5, intStatorMomentumIndex + iCnt + 3);
708  WorkMat.PutColIndex(iCnt + 5, intStatorPositionIndex + iCnt + 3);
709  }
710 
711  WorkMat.PutRowIndex(9, intOmegaRowIndex);
712  WorkMat.PutColIndex(9, intOmegaColumnIndex);
713 
714  Vec3 e3 = m_pStatorNode->GetRCurr().GetCol(3); // corrected axis of rotation of the stator node
715  Vec3 e3_0 = m_pStatorNode->GetRRef().GetCol(3); // predicted axis of rotation of the stator node
716  const Vec3& omega1 = m_pRotorNode->GetWCurr(); // corrected angular velocity of the rotor node
717  const Vec3& omega2 = m_pStatorNode->GetWCurr(); // corrected angular velocity of the rotor node
718  const Vec3& omega1_ref = m_pRotorNode->GetWRef(); // predicted angular velocity of the rotor node
719  const Vec3& omega2_ref = m_pStatorNode->GetWRef(); // predicted angular velocity of the stator node
720 
721  const doublereal OmegaS = m_OmegaS.dGet(); // synchronous angular velocity
722 
723  const doublereal y1 = XCurr(intTorqueDerivativeRowIndex);
724  const doublereal y2 = XCurr(intTorqueRowIndex);
725  const doublereal y5 = XCurr(intOmegaRowIndex);
726 #if 0 // unused
727  const doublereal y1_dot = XPrimeCurr(intTorqueDerivativeRowIndex);
728  const doublereal y2_dot = XPrimeCurr(intTorqueRowIndex);
729 #endif
730  const doublereal y5_dot = XPrimeCurr(intOmegaRowIndex);
731 
732  doublereal s = 1 - y5 / OmegaS; // slip of the rotor
733 
734  if ( std::abs(s) < sm_SingTol )
735  {
736  silent_cerr("\nasynchronous_machine(" << GetLabel() << "): warning slip rate s = " << s << " is close to machine precision!\n");
737  //FIXME: avoid division by zero
738  //FIXME: results might be inaccurate
739  s = copysign(sm_SingTol, s);
740  }
741 
742  doublereal df1_dy1, df1_dy2, df1_dy5;
743  doublereal df1_dy1_dot, df1_dy2_dot, df1_dy5_dot;
744  doublereal df2_dy1, df2_dy2_dot;
745 
746  if (IsMotorOn()) {
747  // power supply is turned on
748 
749  // df1_dy1 = diff(f1,y1)
750  df1_dy1 = ( 2 * m_sK + copysign(1., OmegaS) * y5_dot / ( s * std::pow(OmegaS,2) ) ) * abs(OmegaS);
751  // df1_dy2 = diff(f1,y2)
752  df1_dy2 = ( std::pow(m_sK,2) + std::pow(s,2) ) * std::pow(OmegaS,2) + copysign(1., OmegaS) * y5_dot * m_sK / s;
753  // df1_dy1_dot = diff(f1,diff(y1,t))
754  df1_dy1_dot = 1.;
755  // df2_dy1 = diff(f2,dy1)
756  df2_dy1 = -1.;
757  // df2_dy2_dot = diff(f2,diff(y2,t))
758  df2_dy2_dot = 1.;
759  // df1_ds = diff(f1,s)
760  const doublereal df1_ds = -y1 * y5_dot / ( std::pow(s,2) * OmegaS ) + ( 2 * s * std::pow(OmegaS,2) - copysign(1., OmegaS) * y5_dot * m_sK / std::pow(s,2) ) * y2 - 2 * m_MK * m_sK * std::pow(OmegaS,2);
761  // df1_dy5 = diff(f1,y5) = diff(f1,s) * diff(s,y5)
762  df1_dy5 = df1_ds * ( -1. / OmegaS );
763  // df1_dy5_dot = diff(f1,diff(y5,t))
764  df1_dy5_dot = y1 / ( s * OmegaS ) + copysign(1., OmegaS) * y2 * m_sK / s;
765 
766  } else {
767  // power supply is turned off
768  df1_dy1 = 0.;
769  df1_dy2 = 1.;
770  df1_dy5 = 0.;
771  df1_dy1_dot = 0.;
772  df1_dy2_dot = 0.;
773  df1_dy5_dot = 0.;
774 
775  df2_dy1 = 1.;
776  df2_dy2_dot = 0;
777  }
778 
779  // df3_dy2 = diff(f3,y2)
780  const Vec3 df3_dy2 = e3 * copysign(1., OmegaS); // df3_dy2 = R2 * e3 * sign(OmegaS)
781  // df4_dy2 = diff(f4,y2)
782  const Vec3 df4_dy2 = -df3_dy2; // df4_dy2 = -R2 * e3
783  // df3_dy4 = diff(f3,y4)
784  const Mat3x3 df3_dy4 = -Mat3x3(MatCross, e3_0 * (y2 * copysign(1., OmegaS)) ); // df3_dy4 = -skew( R2^(0) * e3 * y2 * sign(OmegaS) )
785  // df4_dy4 = diff(f4,y4)
786  const Mat3x3 df4_dy4 = -df3_dy4;
787  // df5_dy5 = diff(f5,y5)
788  const doublereal df5_dy5 = 1.;
789  // df5_dy3_dot_T = transpose(diff(f5,diff(y3,t)))
790  const Vec3 df5_dy3_dot_T = -e3; // diff(y5,diff(y3,t)) = -e3^T * R2^T
791  // df5_dy4_dot_T = transpose(diff(f5,diff(y4,t)))
792  const Vec3 df5_dy4_dot_T = -df5_dy3_dot_T; // diff(y5,diff(y4,t)) = e3^T * R2^T
793  // df5_dy3_T = transpose(diff(y5,y3))
794  const Vec3 df5_dy3_T = -omega1_ref.Cross(e3);
795  // df5_dy4_T = transpose(diff(y5,y4))
796  const Vec3 df5_dy4_T = ( omega1 - omega2 ).Cross( e3_0 ) + omega2_ref.Cross( e3 );
797 
798  WorkMat.PutCoef( 1, 1, -df1_dy1_dot - dCoef * df1_dy1 );
799  WorkMat.PutCoef( 1, 2, - dCoef * df1_dy2 );
800  WorkMat.PutCoef( 1, 9, -df1_dy5_dot - dCoef * df1_dy5 );
801  WorkMat.PutCoef( 2, 1, - dCoef * df2_dy1 );
802  WorkMat.PutCoef( 2, 2, -df2_dy2_dot );
803  WorkMat.Put( 3, 2, df3_dy2 * ( -dCoef ) );
804  WorkMat.Put( 3, 6, df3_dy4 * ( -dCoef ) );
805  WorkMat.Put( 6, 2, df4_dy2 * ( -dCoef ) );
806  WorkMat.Put( 6, 6, df4_dy4 * ( -dCoef ) );
807  WorkMat.PutT( 9, 3, -df5_dy3_dot_T + df5_dy3_T * ( -dCoef ) );
808  WorkMat.PutT( 9, 6, -df5_dy4_dot_T + df5_dy4_T * ( -dCoef ) );
809  WorkMat.PutCoef( 9, 9, df5_dy5 * ( -dCoef ) );
810 
811 #ifdef DEBUG
812  std::cerr << __FILE__ ":" << __LINE__ << ":" << __PRETTY_FUNCTION__ << std::endl;
813  std::cerr << "s = " << s << std::endl;
814  std::cerr << "WorkMat=" << std::endl
815  << WorkMat << std::endl << std::endl;
816 #endif
817 
818  return WorkMat_;
819 }
820 
872  doublereal dCoef,
873  const VectorHandler& XCurr,
874  const VectorHandler& XPrimeCurr)
875 {
876 #ifdef DEBUG
877  std::cerr << __FILE__ << ':' << __LINE__ << ':' << __PRETTY_FUNCTION__ << std::endl;
878 #endif
879  integer iNumRows = 0;
880  integer iNumCols = 0;
881 
882  WorkSpaceDim(&iNumRows, &iNumCols);
883 
884  WorkVec.ResizeReset(iNumRows);
885 
886  const integer iFirstIndex = iGetFirstIndex();
887 
888  const integer intTorqueDerivativeRowIndex = iFirstIndex + 1;
889  const integer intTorqueRowIndex = iFirstIndex + 2;
890  const integer intOmegaRowIndex = iFirstIndex + 3;
891 
892  const integer intRotorMomentumIndex = m_pRotorNode->iGetFirstMomentumIndex();
893  const integer intStatorMomentumIndex = m_pStatorNode->iGetFirstMomentumIndex();
894 
895  WorkVec.PutRowIndex(1, intTorqueDerivativeRowIndex);
896  WorkVec.PutRowIndex(2, intTorqueRowIndex);
897 
898  for (int iCnt = 1; iCnt <= 3; ++iCnt) {
899  WorkVec.PutRowIndex(iCnt + 2, intRotorMomentumIndex + iCnt + 3);
900  WorkVec.PutRowIndex(iCnt + 5, intStatorMomentumIndex + iCnt + 3);
901  }
902 
903  WorkVec.PutRowIndex(9, intOmegaRowIndex);
904 
905  Vec3 e3 = m_pStatorNode->GetRCurr().GetCol(3);
906  const Vec3& omega1 = m_pRotorNode->GetWCurr();
907  const Vec3& omega2 = m_pStatorNode->GetWCurr();
908 
909  const doublereal OmegaS = m_OmegaS.dGet();
910 
911  const doublereal y1 = XCurr(intTorqueDerivativeRowIndex); // y1 = diff(M,t)
912  const doublereal y2 = XCurr(intTorqueRowIndex); // y2 = M
913  const doublereal y5 = XCurr(intOmegaRowIndex); // y5 = omega
914  const doublereal y1_dot = XPrimeCurr(intTorqueDerivativeRowIndex); // diff(y1,t) = diff(M,t,2)
915  const doublereal y2_dot = XPrimeCurr(intTorqueRowIndex); // diff(y2,t) = diff(M,t)
916  const doublereal y5_dot = XPrimeCurr(intOmegaRowIndex); // diff(y5,t) = diff(omega,t)
917 
918  doublereal s = 1 - y5 / OmegaS;
919 
920  if ( std::abs(s) < sm_SingTol )
921  {
922  silent_cerr("\nasynchronous_machine(" << GetLabel() << "): warning slip rate s = " << s << " is close to machine precision!\n");
923  //FIXME: avoid division by zero
924  //FIXME: results might be inaccurate
925  s = copysign(sm_SingTol, s);
926  }
927 
928  doublereal f1, f2;
929 
930  if (IsMotorOn()) {
931  // power supply is switched on
932  f1 = y1_dot + ( 2 * m_sK + copysign(1., OmegaS) * y5_dot / ( s * std::pow(OmegaS,2) ) ) * y1 * abs(OmegaS)
933  + ( ( std::pow(m_sK,2) + std::pow(s,2) ) * std::pow(OmegaS,2) + copysign(1., OmegaS) * y5_dot * m_sK / s ) * y2 - 2 * m_MK * m_sK * s * std::pow(OmegaS,2);
934  f2 = y2_dot - y1;
935 
936  } else {
937  // power supply is switched off: torque must be zero
938  f1 = y2;
939  f2 = y1;
940  }
941 
942  const Vec3 f3 = e3 * (y2 * copysign(1., OmegaS));
943  // const Vec3 f4 = -f3;
944 
945  const doublereal f5 = y5 - e3.Dot( omega1 - omega2 );
946 
947  WorkVec.PutCoef( 1, f1 );
948  WorkVec.PutCoef( 2, f2 );
949  WorkVec.Put( 3, f3 );
950  WorkVec.Sub( 6, f3 );
951  WorkVec.PutCoef( 9, f5 );
952 
953 #ifdef DEBUG
954  std::cerr
955  << __FILE__ ":" << __LINE__ << ":" << __PRETTY_FUNCTION__ << std::endl
956  << "s = " << s << std::endl
957  << "y1 = " << y1 << std::endl
958  << "y1_dot = " << y1_dot << std::endl
959  << "y2 = " << y2 << std::endl
960  << "y2_dot = " << y2_dot << std::endl
961  << "f1 = " << f1 << std::endl
962  << "f2 = " << f2 << std::endl
963  << "f3 = " << f3 << std::endl
964  << "f4 = " << -f3 << std::endl
965  << "f5 = " << f5 << std::endl
966  << "WorkVec=" << std::endl
967  << WorkVec << std::endl << std::endl;
968 #endif
969 
970  return WorkVec;
971 }
976 int
978 {
979  return 2;
980 }
981 
986 void
987 asynchronous_machine::GetConnectedNodes(std::vector<const Node *>& connectedNodes) const
988 {
989  connectedNodes.resize(iGetNumConnectedNodes());
990  connectedNodes[0] = m_pRotorNode;
991  connectedNodes[1] = m_pStatorNode;
992 }
993 
1001 void
1003  VectorHandler& X, VectorHandler& XP,
1005 {
1006  const integer intFirstIndex = iGetFirstIndex();
1007 
1008  // 1 2 3
1009  X.Put( intFirstIndex + 1, Vec3(m_dM_dt, m_M, m_omega) );
1010  XP.Put(intFirstIndex + 1, Vec3( 0., m_dM_dt, 0.) );
1011 }
1012 
1019 std::ostream&
1020 asynchronous_machine::Restart(std::ostream& out) const
1021 {
1022  // FIXME: mbdyn crashes when "make restart file;" is specified in the input file before this member function is invoked!
1023  out << "asynchronous_machine, "
1024  "rotor, " << m_pRotorNode->GetLabel() << ", "
1025  "stator, " << m_pStatorNode->GetLabel() << ", "
1026  "MK, " << m_MK << ", "
1027  "sK, " << m_sK << ", "
1028  "OmegaS, ";
1029  m_OmegaS.pGetDriveCaller()->Restart(out) << ", "
1030  "motor on, ";
1031  m_MotorOn.pGetDriveCaller()->Restart(out) << ", "
1032  "M0, " << m_M * copysign(1., m_OmegaS.dGet()) << ", "
1033  "MP0, " << m_dM_dt * copysign(1., m_OmegaS.dGet()) << ";" << std::endl;
1034 
1035  return out;
1036 }
1037 
1045 unsigned int
1047 {
1048  return 0;
1049 }
1050 
1054 void
1056  integer* piNumRows,
1057  integer* piNumCols) const
1058 {
1059  *piNumRows = 0;
1060  *piNumCols = 0;
1061 }
1062 
1068  VariableSubMatrixHandler& WorkMat,
1069  const VectorHandler& XCurr)
1070 {
1071  WorkMat.SetNullMatrix();
1072 
1073  return WorkMat;
1074 }
1075 
1081  SubVectorHandler& WorkVec,
1082  const VectorHandler& XCurr)
1083 {
1084  WorkVec.ResizeReset(0);
1085 
1086  return WorkVec;
1087 }
1088 
1089 bool
1091 {
1092 #ifdef DEBUG
1093  std::cerr << __FILE__ <<":"<< __LINE__ << ":"<< __PRETTY_FUNCTION__ << std::endl;
1094 #endif
1095 
1097 
1098  if (!SetUDE("asynchronous_machine", rf)) {
1099  delete rf;
1100  return false;
1101  }
1102 
1103  return true;
1104 }
1105 
1106 #ifndef STATIC_MODULES
1107 
1108 extern "C"
1109 {
1110 
1116 int
1117 module_init(const char *module_name, void *pdm, void *php)
1118 {
1119  if (!asynchronous_machine_set()) {
1120  silent_cerr("asynchronous_machine: "
1121  "module_init(" << module_name << ") "
1122  "failed" << std::endl);
1123  return -1;
1124  }
1125 
1126  return 0;
1127 }
1128 
1129 } // extern "C"
1130 
1131 #endif // ! STATIC_MODULES
1132 
flag fReadOutput(MBDynParser &HP, const T &t) const
Definition: dataman.h:1064
virtual VariableSubMatrixHandler & InitialAssJac(VariableSubMatrixHandler &WorkMat, const VectorHandler &XCurr)
void PutColIndex(integer iSubCol, integer iCol)
Definition: submat.h:325
virtual unsigned int iGetNumPrivData(void) const
This member function is called when a bind statement or a element drive is used to access private dat...
Vec3 Cross(const Vec3 &v) const
Definition: matvec3.h:218
long int flag
Definition: mbdyn.h:43
virtual const Mat3x3 & GetRRef(void) const
Definition: strnode.h:1006
virtual bool bToBeOutput(void) const
Definition: output.cc:890
GradientExpression< BinaryExpr< FuncPow, LhsExpr, RhsExpr > > pow(const GradientExpression< LhsExpr > &u, const GradientExpression< RhsExpr > &v)
Definition: gradient.h:2961
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
Definition: matvec3.h:98
virtual void ResizeReset(integer)
Definition: vh.cc:55
const MatCross_Manip MatCross
Definition: matvec3.cc:639
FullSubMatrixHandler & SetFull(void)
Definition: submat.h:1168
virtual const Mat3x3 & GetRCurr(void) const
Definition: strnode.h:1012
void PutT(integer iRow, integer iCol, const Vec3 &v)
Definition: submat.cc:239
asynchronous_machine(unsigned uLabel, const DofOwner *pDO, DataManager *pDM, MBDynParser &HP)
Constructs the element from the information in the input file.
bool IsMotorOn(void) const
If a drive caller has been specified in the input file the value returned by the drive caller is used...
doublereal Dot(const Vec3 &v) const
Definition: matvec3.h:243
SubVectorHandler & AssRes(SubVectorHandler &WorkVec, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
This member function implements the equation of an asynchronous machine according to Maschinendynamik...
virtual void Sub(integer iRow, const Vec3 &v)
Definition: vh.cc:78
void Put(integer iRow, integer iCol, const Vec3 &v)
Definition: submat.cc:221
void PutCoef(integer iRow, integer iCol, const doublereal &dCoef)
Definition: submat.h:672
Vec3 GetCol(unsigned short int i) const
Definition: matvec3.h:903
VariableSubMatrixHandler & AssJac(VariableSubMatrixHandler &WorkMat, doublereal dCoef, const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
Computes the contribution to the jacobian matrix .
virtual DofOrder::Order GetDofType(unsigned int i) const
The type of the equation of the private degrees of freedom is determined by this member function...
virtual const Vec3 & GetWRef(void) const
Definition: strnode.h:1024
std::vector< Hint * > Hints
Definition: simentity.h:89
virtual void Update(const VectorHandler &XCurr, const VectorHandler &XPrimeCurr)
This member function is called after each iteration in the nonlinear solution phase.
virtual std::ostream & Restart(std::ostream &out) const =0
virtual void Output(OutputHandler &OH) const
Writes private data to a file at each time step.
virtual void PutRowIndex(integer iSubRow, integer iRow)=0
virtual int iGetNumConnectedNodes(void) const
This member function is called if the statement "print: node connection;" is specified in the input f...
virtual doublereal dGetPrivData(unsigned int i) const
This member function is called whenever a bind statement or a element drive is used to access private...
int module_init(const char *module_name, void *pdm, void *php)
This function registers our user defined element for the math parser.
void SetNullMatrix(void)
Definition: submat.h:1159
virtual bool IsKeyWord(const char *sKeyWord)
Definition: parser.cc:910
bool asynchronous_machine_set(void)
std::ostream & Loadable(void) const
Definition: output.h:506
doublereal copysign(doublereal x, doublereal y)
Definition: gradient.h:97
virtual SubVectorHandler & InitialAssRes(SubVectorHandler &WorkVec, const VectorHandler &XCurr)
virtual std::ostream & DescribeDof(std::ostream &out, const char *prefix, bool bInitial) const
This member function is called if the statement "print: dof description;" is specified in the input f...
virtual integer iGetFirstMomentumIndex(void) const =0
virtual integer iGetFirstPositionIndex(void) const
Definition: strnode.h:452
virtual const Vec3 & GetWCurr(void) const
Definition: strnode.h:1030
unsigned int uLabel
Definition: withlab.h:44
virtual void SetInitialValue(VectorHandler &X)
virtual void PutCoef(integer iRow, const doublereal &dCoef)=0
VectorExpression< VectorCrossExpr< VectorLhsExpr, VectorRhsExpr >, 3 > Cross(const VectorExpression< VectorLhsExpr, 3 > &u, const VectorExpression< VectorRhsExpr, 3 > &v)
Definition: matvec.h:3248
DriveCaller * pGetDriveCaller(void) const
Definition: drive.cc:658
virtual void GetConnectedNodes(std::vector< const Node * > &connectedNodes) const
This member function is called if the statement "print: node connection;" is specified in the input f...
virtual void ResizeReset(integer, integer)
Definition: submat.cc:182
Definition: except.h:79
virtual void AfterPredict(VectorHandler &X, VectorHandler &XP)
virtual void Put(integer iRow, const Vec3 &v)
Definition: vh.cc:93
virtual bool IsArg(void)
Definition: parser.cc:807
Definition: elem.h:75
static const doublereal sm_SingTol
void PutRowIndex(integer iSubRow, integer iRow)
Definition: submat.h:311
std::ostream & GetLogFile(void) const
Definition: dataman.h:326
bool SetUDE(const std::string &s, UserDefinedElemRead *rude)
Definition: userelem.cc:97
void Set(const DriveCaller *pDC)
Definition: drive.cc:647
doublereal dGet(const doublereal &dVar) const
Definition: drive.cc:664
virtual void WorkSpaceDim(integer *piNumRows, integer *piNumCols) const
The size of the contribution to the residual and the jacobian matrix is determined by this member fun...
virtual std::ostream & DescribeEq(std::ostream &out, const char *prefix, bool bInitial) const
This member function is called if the statement "print: equation description;" is specified in the in...
DriveCaller * GetDriveCaller(bool bDeferred=false)
Definition: mbpar.cc:2033
virtual unsigned int iGetInitialNumDof(void) const
This member function must be implemented only if the initial assembly feature is requested.
virtual std::ostream & Restart(std::ostream &out) const
This member function is called if "make restart file;" is specified in the input file.
virtual void SetOutputFlag(flag f=flag(1))
Definition: output.cc:896
virtual integer iGetFirstIndex(void) const
Definition: dofown.h:127
virtual unsigned int iGetPrivDataIdx(const char *s) const
The following string values can be specified in a bind statement or in an element drive caller...
double doublereal
Definition: colamd.c:52
long int integer
Definition: colamd.c:51
virtual HighParser::ErrOut GetLineData(void) const
Definition: parsinc.cc:697
unsigned int GetLabel(void) const
Definition: withlab.cc:62
Node * ReadNode(MBDynParser &HP, Node::Type type) const
Definition: dataman3.cc:2309
virtual unsigned int iGetNumDof(void) const
The number of private degrees of freedom is determined by this member function.
bool UseText(int out) const
Definition: output.cc:446
virtual void SetValue(DataManager *pDM, VectorHandler &X, VectorHandler &XP, SimulationEntity::Hints *ph)
This member function is called before the integration starts in order to set the initial values for t...
virtual void InitialWorkSpaceDim(integer *piNumRows, integer *piNumCols) const
virtual doublereal GetReal(const doublereal &dDefval=0.0)
Definition: parser.cc:1056