added SIM_singlerigidbody (#1762)

* added SIM_singlerigidbody

* Removed a few sections in the README file

* Removed the sims directory

* Added shutdown job and made minor changes to body.cpp

---------

Co-authored-by: Wallace <bnwalla1@scooby.trick.gov>
This commit is contained in:
bnwalla1 2024-10-16 14:37:48 -05:00 committed by GitHub
parent b2403dcfc6
commit 7bc6ab74a8
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
27 changed files with 709 additions and 0 deletions

Binary file not shown.

After

Width:  |  Height:  |  Size: 1.3 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 3.1 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 3.4 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 1.8 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 6.9 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 1.6 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 886 B

Binary file not shown.

After

Width:  |  Height:  |  Size: 3.3 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 1.8 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 2.0 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 2.9 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 3.4 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 1.7 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 2.8 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 1.5 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 3.2 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 157 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 88 KiB

View File

@ -0,0 +1,12 @@
trick.frame_log_on()
trick.real_time_enable()
trick.exec_set_software_frame(0.033)
trick.itimer_enable()
trick.exec_set_enable_freeze(True)
#trick.exec_set_freeze_command(True)
#simControlPanel = trick.SimControlPanel()
#trick.add_external_application(simControlPanel)

View File

@ -0,0 +1,54 @@
global DR_GROUP_ID
global drg
try:
if DR_GROUP_ID >= 0:
DR_GROUP_ID += 1
except NameError:
DR_GROUP_ID = 0
drg = []
drg.append(trick.DRBinary("StateData"))
drg[DR_GROUP_ID].set_freq(trick.DR_Always)
drg[DR_GROUP_ID].set_cycle(0.10)
drg[DR_GROUP_ID].set_single_prec_only(False)
drg[DR_GROUP_ID].add_variable("dyn.body.force[0]")
drg[DR_GROUP_ID].add_variable("dyn.body.force[1]")
drg[DR_GROUP_ID].add_variable("dyn.body.force[2]")
drg[DR_GROUP_ID].add_variable("dyn.body.angle_force[0]")
drg[DR_GROUP_ID].add_variable("dyn.body.angle_force[1]")
drg[DR_GROUP_ID].add_variable("dyn.body.angle_force[2]")
drg[DR_GROUP_ID].add_variable("dyn.body.rotate[0][0]")
drg[DR_GROUP_ID].add_variable("dyn.body.rotate[0][1]")
drg[DR_GROUP_ID].add_variable("dyn.body.rotate[0][2]")
drg[DR_GROUP_ID].add_variable("dyn.body.rotate[1][0]")
drg[DR_GROUP_ID].add_variable("dyn.body.rotate[1][1]")
drg[DR_GROUP_ID].add_variable("dyn.body.rotate[1][2]")
drg[DR_GROUP_ID].add_variable("dyn.body.rotate[2][0]")
drg[DR_GROUP_ID].add_variable("dyn.body.rotate[2][1]")
drg[DR_GROUP_ID].add_variable("dyn.body.rotate[2][2]")
drg[DR_GROUP_ID].add_variable("dyn.body.pos[0]")
drg[DR_GROUP_ID].add_variable("dyn.body.pos[1]")
drg[DR_GROUP_ID].add_variable("dyn.body.pos[2]")
drg[DR_GROUP_ID].add_variable("dyn.body.angle[0]")
drg[DR_GROUP_ID].add_variable("dyn.body.angle[1]")
drg[DR_GROUP_ID].add_variable("dyn.body.angle[2]")
drg[DR_GROUP_ID].add_variable("dyn.body.vel[0]")
drg[DR_GROUP_ID].add_variable("dyn.body.vel[1]")
drg[DR_GROUP_ID].add_variable("dyn.body.vel[2]")
drg[DR_GROUP_ID].add_variable("dyn.body.acc[0]")
drg[DR_GROUP_ID].add_variable("dyn.body.acc[1]")
drg[DR_GROUP_ID].add_variable("dyn.body.acc[2]")
drg[DR_GROUP_ID].add_variable("dyn.body.omega[0]")
drg[DR_GROUP_ID].add_variable("dyn.body.omega[1]")
drg[DR_GROUP_ID].add_variable("dyn.body.omega[2]")
drg[DR_GROUP_ID].add_variable("dyn.body.omegaDot[0]")
drg[DR_GROUP_ID].add_variable("dyn.body.omegaDot[1]")
drg[DR_GROUP_ID].add_variable("dyn.body.omegaDot[2]")
drg[DR_GROUP_ID].add_variable("dyn.body.FORCE_INIT[0]")
drg[DR_GROUP_ID].add_variable("dyn.body.FORCE_INIT[1]")
drg[DR_GROUP_ID].add_variable("dyn.body.FORCE_INIT[2]")
trick.add_data_record_group(drg[DR_GROUP_ID], trick.DR_Buffer)
drg[DR_GROUP_ID].enable()

View File

@ -0,0 +1,122 @@
# Single Rigid Body
``SIM_singlerigidbody`` is a simulation of the motion of a sphere that is a single rigid body with six degrees of freedom.
The motion of the sphere is found by using Euler's equations of motion, and can be controlled by changing several initial input values such as force, position, angle, angular velocity, velocity, and position of the force. The recorded data of the motion can then be viewed and analyzed.
## Dynamics model
### Equations of Motion
To determine the sphere's motion, we will use Euler's equations of motion for a single rigid body. Euler's equations describe the combined translational and rotational dynamics of a rigid body in three-dimensional space.
![Equations_of_motion](images\EquationsofMotion.png)
Where,
![Force_and_Torque_eq](images\ForceTorque.png)
* is a 6x1 column vector.
* F is the applied force.
* &#0964; is the torque.
![Mass_Matrix_eq](images\MassMatrix.png)
* is a 6x6 matrix.
* m is the mass.
* I <sub> 3x3 </sub> is the identity matrix.
* I is the inertia.
* C <sub> m </sub> is the center of mass.
![Acceleration_Omega](images\AccOmega.png)
* is a 6x1 column vector.
* a<sub> o </sub> is the acceleration.
![Force_Torque_Remain](images\ForceTorqueRemain.png)
* is a 6x1 column vector.
* &#0969; <sup> A </sup> is the angular velocity.
#### Force
The force acting on the sphere is calculated by the cross product of the rotation matrix and initial force.
![Force](images\Force.png)
#### Torque
Torque is found by the cross product of the rotation matrix and the cross product of the position of force and initial force.
![Torque](images\Torque.png)
#### Mass
The default value of the mass of the sphere is:
![Mass](images\Mass.png)
#### Radius
The default value of the radius is:
![Radius](images\Radius.png)
#### Center of mass
The center of mass of a sphere is equal to its radius.
![Center_of_mass](images\Centerofmass.png)
#### Inertia
Inertia of a sphere is calculated by:
![Inertia](images\Inertia.png)
#### Acceleration and angular acceleration
To determine the acceleration and angular acceleration, ```dLU_Choleski``` must be used. ```dLU_Choleski``` solves linear sets of equations in the form of [A]x = b and returns the vector x, which is the acceleration. It decomposes the symmetric, positive definite [A] matrix into a lower triangular form, where we can then substitute and our equation becomes:
![Vector_a](images\Vectora.png)
The term 'a' is a 6x1 vector. Acceleration is equal to the first three values of vector a and angular acceleration is equal to the last three values of vector a.
### State Integration
The state is defined by position, velocity, angle, and angular velocity. They are calculated by numerically integrating velocity, acceleration, angular velocity, and angular acceleration, respectively.
#### Position
The position of the sphere is determined by integrating the velocity over time.
![Position](images\Position.png)
#### Velocity
The velocity is determined by integrating the acceleration of the sphere over time.
![Velocity](images\Velocity.png)
#### Angle
The angle is determined by integrating the angular velocity over time.
![Angle](images\Angle.png)
#### Angular velocity
The angular velocity is determined by integrating the angular acceleration over time.
![Angular_Velocity](images\Angularvelocity.png)
## Inputs
Variable |Type |Units
----------------------------|-------------|---------
dyn.body.POS_INIT |double[3] |m
dyn.body.ANGLE_INIT |double[3] |rad
dyn.body.OMEGA_INIT |double[3] |rad*s^-1
dyn.body.FORCE_INIT |double[3] |N
dyn.body.TORQUE_INIT |double[3] |N*m
dyn.body.VEL_INIT |double[3] |m/s
dyn.body.POS_FORCE |double[3] |m
## Outputs
Variable |Type |Units
----------------------------|-------------|--------
dyn.body.pos |double[3] |m
dyn.body.vel |double[3] |m/s
dyn.body.acc |double[3] |m/s^2
dyn.FORCE_INIT |double[3] |N
dyn.force |double[3] |N
dyn.body.angle_force |double[3] |rad
dyn.body.omegaDot |double[3] |rad/^s2
dyn.body.omega |double[3] |rad*s^-1
dyn.body.rotate |double[3] |

View File

@ -0,0 +1,81 @@
################################################################################
# _ ___ _ #
# / |___| _ ) ___ __| |_ _ #
# | |___| _ \/ _ \/ _` | || | #
# |_| |___/\___/\__,_|\_, | #
# __ ___ ___ __|__/ #
# / / ___| \ / _ \| __| #
# / _ \___| |) | (_) | _| #
# \___/ |___/ \___/|_| #
# #
# This input creates a simple body with 6 DOF: #
# -Creates 1 body #
# -Initializes 6 DOF #
################################################################################
# ============================================================================ -
# This file is used to store several quality of life python functions.
# ============================================================================ -
exec(open("Modified_data/state_data.dr").read())
exec(open("Modified_data/realtime.py").read())
################################################################################
# dyn DATA #
# (NOTE: dyn is our instance of BodySimObject) #
# This builds our objects, and defines their physical properties. #
################################################################################
# ============================================================================ -
# Establish the physical properties of the body...CONSTANTS!
# ============================================================================ -
dyn.body.mass = 20.0
dyn.body.radius = 2.0
# ============================================================================ -
# Initialize system
# ============================================================================ -
M_PI = 3.141592653589793238460
#dyn.body.POS_INIT[0] = 0.0
#dyn.body.POS_INIT[1] = 0.0
#dyn.body.POS_INIT[2] = 0.0
#dyn.body.ANGLE_INIT[0] = 0.0
#dyn.body.ANGLE_INIT[1] = 0.0
#dyn.body.ANGLE_INIT[2] = 0.0
#dyn.body.OMEGA_INIT[0] = 0.0
#dyn.body.OMEGA_INIT[1] = 0.0
#dyn.body.OMEGA_INIT[2] = 0.0
#dyn.body.FORCE_INIT[0] = 0.0
#dyn.body.FORCE_INIT[1] = 0.0
#dyn.body.FORCE_INIT[2] = 10.0
#dyn.body.TORQUE_INIT[0] = 0.0
#dyn.body.TORQUE_INIT[1] = 0.0
#dyn.body.TORQUE_INIT[2] = 0.0
#dyn.body.VEL_INIT[0] = 0.0
#dyn.body.VEL_INIT[1] = 0.0
#dyn.body.VEL_INIT[2] = 0.0
#dyn.body.POS_FORCE[0] = 2.0
#dyn.body.POS_FORCE[1] = 0.0
#dyn.body.POS_FORCE[2] = 0.0
################################################################################
# SIMULATION SET UP #
# This is where we set up the remaining inputs to create a working sim. #
################################################################################
# ============================================================================ -
# getintegrator(integration method, number of variables to be integrated)
# This sets up the integration scheme we use to analyze this system.
# ============================================================================ -
dyn_integloop.getIntegrator(trick.Runge_Kutta_4, 12)
# ============================================================================ -
# This tells the sim to stop after a few seconds
# ============================================================================ -
trick.stop(10.0)

View File

@ -0,0 +1,32 @@
/************************TRICK HEADER*************************
PURPOSE:
( Single body rigid dynamics simulation )
LIBRARY DEPENDENCIES:
((src/body.cpp)
(src/body_shutdown.cpp)
)
*************************************************************/
#define TRICK_NO_MONTE_CARLO
#define TRICK_NO_MASTERSLAVE
#define TRICK_NO_INSTRUMENTATION
#define TRICK_NO_REALTIMEINJECTOR
#define TRICK_NO_ZEROCONF
#include "sim_objects/default_trick_sys.sm"
##include "include/body.hh"
class BodySimObject : public Trick::SimObject {
public:
BODY body;
BodySimObject() {
("default_data") body.default_data() ;
("initialization") body.init() ;
("derivative") body.derivative() ;
("integration") trick_ret = body.integ() ;
("shutdown") body.body_shutdown();
}
};
BodySimObject dyn;
IntegLoop dyn_integloop(0.01) dyn;

View File

@ -0,0 +1,2 @@
TRICK_CFLAGS += -Imodels
TRICK_CXXFLAGS += -Imodels

View File

@ -0,0 +1,78 @@
/*********************************************************************
PURPOSE: (Simulate a single rigid body with six degrees of freedom)
LIBRARY DEPENDENCIES:
((SingleRigidBody/src/body.o))
*********************************************************************/
#ifndef BODY_HH
#define BODY_HH
class BODY {
public:
BODY();
~BODY();
//Variables
//Initializes state variables
double POS_INIT[3];
double VEL_INIT[3];
double ANGLE_INIT[3];
double OMEGA_INIT[3];
//State variables
double pos[3];
double vel[3];
double angle[3];
double omega[3];
double acc[3];
double omegaDot[3];
double force[3];
double FORCE_INIT[3];
double force_remain[3];
double POS_FORCE[3];
double torque[3];
double TORQUE_INIT[3];
double torque_remain[3];
double inertia;
double inertia_matrix[3][3];
double radius;
double mass;
double mat_mass[6][6];
double massmatrix[3][3];
double CM[3];
double CM_skew[3][3];
double m_CM_skew[3][3];
double neg_m_CM_skew[3][3];
double vec_a[6];
double vec_b[6];
double rotate[3][3];
double angle_force[3];
double tmp_vec[6];
double **mat_mass_dyn;
double **mat_L;
// Methods
void default_data();
void init();
void derivative();
int integ();
int body_shutdown();
void rotation_matrix();
void calcforce();
void calctorque();
void calcforce_remain();
void calctorque_remain();
void eq_setup();
void eq_solver();
};
#endif

View File

@ -0,0 +1,314 @@
/**************************************************************************
PURPOSE: (Test trick macros (matrix and vector) as well as trick functions
including euler_matrix, deuler_123, and dLU_Choleski. Perform a
single body motion with six degrees of freedom)
LIBRARY DEPENDENCIES:
((src/body.o))
**************************************************************************/
#include <math.h>
#include "trick/integrator_c_intf.h"
#include "trick/trick_math_proto.h"
#include "trick/trick_math_error.h"
#include "trick/reference_frame.h"
#include "trick/vector_macros.h"
#include "trick/matrix_macros.h"
#include "../include/body.hh"
BODY::BODY(){};
BODY::~BODY(){};
/*************************************************************************
FUNCTION: SingleRigidBody::default_data()
PURPOSE:` (Creates an object)
*************************************************************************/
void BODY::default_data() {
//Initializing vectors and matricies to zero.
for(int i = 0; i<6; i++)
vec_a[i] = 0.0;
for(int i = 0; i<6; i++)
vec_b[i] = 0.0;
for(int i = 0; i<6; i++)
for(int j = 0; j<6; j++)
mat_mass[i][j] = 0.0;
V_INIT(pos);
V_INIT(vel);
V_INIT(acc);
V_INIT(force);
V_INIT(force_remain);
V_INIT(torque);
V_INIT(torque_remain);
V_INIT(angle);
V_INIT(omega);
V_INIT(omegaDot);
V_INIT(angle_force);
M_INIT(CM_skew);
M_INIT(m_CM_skew);
M_INIT(massmatrix);
M_INIT(rotate);
// Inputs
mass = 20.0;
radius = 2.0;
POS_INIT[0] = 0.0;
POS_INIT[1] = 0.0;
POS_INIT[2] = 0.0;
ANGLE_INIT[0] = 0.0;
ANGLE_INIT[1] = 0.0;
ANGLE_INIT[2] = 0.0;
OMEGA_INIT[0] = 0.0;
OMEGA_INIT[1] = 0.0;
OMEGA_INIT[2] = 0.0;
FORCE_INIT[0] = 0.0;
FORCE_INIT[1] = 0.0;
FORCE_INIT[2] = 10.0;
VEL_INIT[0] = 0.0;
VEL_INIT[1] = 0.0;
VEL_INIT[2] = 0.0;
POS_FORCE[0] = 2.0;
POS_FORCE[1] = 0.0;
POS_FORCE[2] = 0.0;
};
/************************************************************************
FUNCTION: SingleRigidBody::init()
PURPOSE: (Initializing state variables, creating inertia and CoM skew
matrices)
************************************************************************/
void BODY::init() {
V_COPY(pos, POS_INIT);
V_COPY(vel, VEL_INIT);
V_COPY(angle, ANGLE_INIT);
V_COPY(omega, OMEGA_INIT);
// Center of mass skew matrix
V_INIT(CM);
V_SKEW(CM_skew, CM);
MxSCALAR(m_CM_skew, CM_skew, mass);
MxSCALAR(neg_m_CM_skew, m_CM_skew, -1.0);
// Creates 3x3 matrix for inerta
inertia = (2.0/5.0) * mass * radius * radius;
inertia_matrix[0][0] = inertia_matrix[1][1] = inertia_matrix[2][2] = inertia;
// Creates 3x3 matrix for mass
massmatrix[0][0] = massmatrix[1][1] = massmatrix[2][2] = mass;
// Combines inertia, center of mass, and mass matrix into a 6x6 matrix
for(int i = 0; i<3; i++)
for(int j = 0; j<3; j++)
mat_mass[i][j] = massmatrix[i][j];
for(int i = 0; i<3; i++)
for(int j = 3; j<6; j++)
mat_mass[i][j] = neg_m_CM_skew[i][j-3];
for(int i = 3; i<6; i++)
for(int j = 0; j<3; j++)
mat_mass[i][j] = m_CM_skew[i-3][j];
for(int i = 3; i<6; i++)
for(int j = 3; j<6; j++)
mat_mass[i][j] = inertia_matrix[i-3][j-3];
// Temporary vector and dynamic memory for Choleski
for(int i = 0; i<6; i++)
tmp_vec[i] = 0.0;
mat_L = new double*[6];
mat_mass_dyn = new double*[6];
};
/*****************************************************************************
FUNCTION: BODY::rotation_matrix()
PURPOSE: (Calculates transformation matrix from provided angles)
*****************************************************************************/
void BODY::rotation_matrix() {
// euler_matrix outputs the rotate matrix from angle
euler_matrix(angle, rotate, 0, Roll_Pitch_Yaw);
};
/******************************************************************************
FUNCTION: BODY::calcforce()
PURPOSE: (Calculates force wrt inertial frame)
******************************************************************************/
void BODY::calcforce() {
MxV(force, rotate, FORCE_INIT);
};
/*****************************************************************************
FUNCTION: BODY::calctorque()
PURPOSE: (Calculates torque wrt inertial frame)
*****************************************************************************/
void BODY::calctorque() {
double torque_force[3];
V_CROSS(torque_force, POS_FORCE, FORCE_INIT);
MxV(torque, rotate, torque_force);
};
/****************************************************************************
FUNCTION: BODY::calcforce_remain()
PURPOSE: (Calculate remaining force wrt inertial frame)
****************************************************************************/
void BODY::calcforce_remain() {
double cross1[3];
double cross2[3];
V_CROSS(cross1, omega, CM);
V_CROSS(cross2, omega, cross1);
V_SCALE(force_remain, cross2, mass);
};
/***************************************************************************
FUNCTION: BODY::calctorque_remain()
PURPOSE: (Calculate remaining torque wrt inertial frame)
***************************************************************************/
void BODY::calctorque_remain() {
double I_w[3];
MxV(I_w, inertia_matrix, omega);
V_CROSS(torque_remain, omega, I_w);
};
/***************************************************************************
FUNCTION: BODY::eq_setup()
PURPOSE: (Establishes b column vector for final eq M*a=b)
***************************************************************************/
void BODY::eq_setup() {
double remainder_force[3];
double remainder_torque[3];
// Subtracts remaining from force & torque to get vec_b.
V_SUB(remainder_force, force, force_remain);
V_SUB(remainder_torque, torque, torque_remain);
for(int i = 0; i<3; i++){
vec_b[i] = remainder_force[i];
}
for(int i = 3; i<6; i++){
vec_b[i] = remainder_torque[i-3];
}
};
/**************************************************************************
FUNCTION: BODY::eq_solver()
PURPOSE: (Performs choleski decomposition and computes linear and angular
accelerations)
**************************************************************************/
void BODY::eq_solver() {
// Solving a = b * M^-1
for(int i = 0; i<6; i++)
{
mat_mass_dyn[i] = mat_mass[i];
mat_L[i] = new double[6];
}
dLU_Choleski(mat_mass_dyn, mat_L, tmp_vec, 6, vec_b, vec_a, 0);
// Acceleration output
for(int i = 0; i<3; i++)
acc[i] = vec_a[i];
// Angular accleration output
for(int i = 0; i<3; i++)
omegaDot[i] = vec_a[i + 3];
};
/***************************************************************************
FUNCTION: BODY::derivative()
PURPOSE: (Calls all functions in desired order for calculations)
***************************************************************************/
void BODY::derivative() {
rotation_matrix();
calcforce();
calctorque();
calcforce_remain();
calctorque_remain();
eq_setup();
eq_solver();
};
/***************************************************************************
FUNCTION: BODY::integrate()
PURPOSE: (Sets up trick integration)
***************************************************************************/
int BODY::integ() {
int integration_step;
load_state(
&pos[0], &pos[1], &pos[2],
&vel[0], &vel[1], &vel[2],
&angle[0], &angle[1], &angle[2],
&omega[0], &omega[1], &omega[2],
NULL);
load_deriv(
&vel[0], &vel[1], &vel[2],
&acc[0], &acc[1], &acc[2],
&omega[0], &omega[1], &omega[2],
&omegaDot[0], &omegaDot[1], &omegaDot[2],
NULL);
integration_step = integrate();
unload_state(
&pos[0], &pos[1], &pos[2],
&vel[0], &vel[1], &vel[2],
&angle[0], &angle[1], &angle[2],
&omega[0], &omega[1], &omega[2],
NULL);
return(integration_step);
};

View File

@ -0,0 +1,14 @@
/************************************************************************
PURPOSE: (Shutdown the simulation)
*************************************************************************/
#include "../include/body.hh"
#include "trick/exec_proto.h"
int BODY::body_shutdown() {
delete[] mat_mass_dyn;
delete[] mat_L;
return 0;
};