diff --git a/src/demos/core/CMakeLists.txt b/src/demos/core/CMakeLists.txt index 7163f4959d..ee7b4238c7 100644 --- a/src/demos/core/CMakeLists.txt +++ b/src/demos/core/CMakeLists.txt @@ -1,7 +1,7 @@ #-------------------------------------------------------------- # Add executables -SET(DEMOS +set(DEMOS demo_CH_coords demo_CH_linalg demo_CH_matrix_ref @@ -15,23 +15,38 @@ SET(DEMOS demo_CH_filesystem ) +#-------------------------------------------------------------- + +# Set include paths, compiler & linker flags, and libraries + +include_directories(${CMAKE_CURRENT_SOURCE_DIR}) + +include_directories(${CH_INCLUDES}) +set(COMPILER_FLAGS "${CH_CXX_FLAGS}") +set(LINKER_FLAGS "${CH_LINKERFLAG_EXE}") +list(APPEND LIBS "ChronoEngine") + +if(ENABLE_MODULE_PARDISO_MKL) + include_directories(${CH_MKL_INCLUDES}) + set(COMPILER_FLAGS "${COMPILER_FLAGS} ${CH_MKL_CXX_FLAGS}") + set(LINKER_FLAGS "${LINKER_FLAGS} ${CH_MKL_LINK_FLAGS}") + list(APPEND LIBS "ChronoEngine_pardisomkl") +endif() + +#-------------------------------------------------------------- -MESSAGE(STATUS "Demo programs for CORE module...") +message(STATUS "Demo programs for CORE module...") -FOREACH(PROGRAM ${DEMOS}) - MESSAGE(STATUS "...add ${PROGRAM}") +foreach(PROGRAM ${DEMOS}) + message(STATUS "...add ${PROGRAM}") - ADD_EXECUTABLE(${PROGRAM} "${PROGRAM}.cpp") - SOURCE_GROUP("" FILES "${PROGRAM}.cpp") + add_executable(${PROGRAM} "${PROGRAM}.cpp") + source_group("" FILES "${PROGRAM}.cpp") - SET_TARGET_PROPERTIES(${PROGRAM} PROPERTIES - FOLDER demos - COMPILE_FLAGS "${CH_CXX_FLAGS}" - LINK_FLAGS "${CH_LINKERFLAG_EXE}") - SET_PROPERTY(TARGET ${PROGRAM} PROPERTY VS_DEBUGGER_WORKING_DIRECTORY "$") - TARGET_LINK_LIBRARIES(${PROGRAM} ChronoEngine) - ADD_DEPENDENCIES(${PROGRAM} ChronoEngine) + set_target_properties(${PROGRAM} PROPERTIES FOLDER demos COMPILE_FLAGS "${COMPILER_FLAGS}" LINK_FLAGS "${LINKER_FLAGS}") + set_property(TARGET ${PROGRAM} PROPERTY VS_DEBUGGER_WORKING_DIRECTORY "$") + target_link_libraries(${PROGRAM} ${LIBS}) - INSTALL(TARGETS ${PROGRAM} DESTINATION ${CH_INSTALL_DEMO}) -ENDFOREACH(PROGRAM) + install(TARGETS ${PROGRAM} DESTINATION ${CH_INSTALL_DEMO}) +endforeach(PROGRAM) diff --git a/src/demos/core/demo_CH_powertrain.cpp b/src/demos/core/demo_CH_powertrain.cpp index bc3edc942d..3c5f83bd3a 100644 --- a/src/demos/core/demo_CH_powertrain.cpp +++ b/src/demos/core/demo_CH_powertrain.cpp @@ -9,7 +9,7 @@ // http://projectchrono.org/license-chrono.txt. // // ============================================================================= -// Authors: Alessandro Tasora +// Authors: Alessandro Tasora, Radu Serban // ============================================================================= // // Demo code about creating a power train using the '1D' items of ChShaft class @@ -34,562 +34,683 @@ #include "chrono_thirdparty/filesystem/path.h" -using namespace chrono; +#include "demos/SetChronoSolver.h" -int main(int argc, char* argv[]) { - std::cout << "Copyright (c) 2017 projectchrono.org\nChrono version: " << CHRONO_VERSION << std::endl; +using namespace chrono; - // Create output directory - std::string out_dir = GetChronoOutputPath() + "DEMO_POWERTRAIN"; - if (!filesystem::create_directory(filesystem::path(out_dir))) { - std::cout << "Error creating directory " << out_dir << std::endl; - return 1; +void Example1(const std::string& out_dir, ChSolver::Type solver_type, ChTimestepper::Type integrator_type) { + std::cout << "\n-------------------------------------------------------------------------\n"; + std::cout << "Example: create a simple power train with ChShaft objects\n" << std::endl; + + // We will model a very basic powertrain with two shafts A and B, + // connected by a reducer [ t ] with transmission ratio 't'. Shafts are + // free to rotate, shaft A has an applied torque Ta, so A and B will + // constantly accelerate. Each shaft must have some inertia, it's like a + // flywheel, marked as || in the following scheme: + // + // A B + // Ta ||---[ t ]---|| + // + + // The physical system: it contains all physical objects. + ChSystemNSC sys; + SetChronoSolver(sys, solver_type, integrator_type); + + // Create a 1-degree-of-freedom '1D' mechanical object, that + // is a ChShaft (an item that can oly rotate, with one inertia value + // and maybe one applied torque). The ChShaft objects do not have + // any meaning in 3d: they are just 'building blocks' for making + // power trains as in imput-output black box schemes. + auto my_shaftA = chrono_types::make_shared(); + my_shaftA->SetInertia(10); + my_shaftA->SetAppliedLoad(6); + sys.AddShaft(my_shaftA); + + // Create another shaft. Note that we use shared pointers for ChShaft + // objects, as we did for ChBody objects. + auto my_shaftB = chrono_types::make_shared(); + my_shaftB->SetInertia(100); + my_shaftB->SetFixed(false); + sys.AddShaft(my_shaftB); + + // Create a ChShaftsGear, that represents a simplified model + // of a reducer, with transmission ratio t, between two ChShaft objects. + // (Note that you could also build a 3D powertrain by creating full rigid bodies + // of ChBody type and join them using ChLinkLockRevolute, ChLinkLockGear 3D constraints, + // but this would introduce many unnecessary degrees of freedom/constraints + // whereas the 1D items of ChShaft type, in this example, make things much simplier). + auto my_shaft_gearAB = chrono_types::make_shared(); + my_shaft_gearAB->Initialize(my_shaftA, my_shaftB); + my_shaft_gearAB->SetTransmissionRatio(-0.1); // ex., a couple of spur gears with 20 and 200 teeth + sys.Add(my_shaft_gearAB); + + std::cout << "\nSystem hierarchy" << std::endl; + sys.ShowHierarchy(std::cout); + + std::ofstream file_results(out_dir + "/example1.txt"); + + // Simulation loop + double end_time = 2.5; + double step_size = 0.01; + + for (double frame_time = 0.05; frame_time < end_time; frame_time += 0.05) { + // Perform simulation up to frame_time + sys.DoFrameDynamics(frame_time, step_size); + + // Print results + std::cout << "Time: " << frame_time << std::endl + << " shaft A rot: " << my_shaftA->GetPos() << " speed: " << my_shaftA->GetPosDt() + << " accel: " << my_shaftA->GetPosDt2() << std::endl + << " shaft B rot: " << my_shaftB->GetPos() << " speed: " << my_shaftB->GetPosDt() + << " accel: " << my_shaftB->GetPosDt2() << std::endl + << " AB gear, torque on A side: " << my_shaft_gearAB->GetReaction1() + << " AB gear, torque on B side: " << my_shaft_gearAB->GetReaction2() << std::endl; + file_results << frame_time << ", " << my_shaftA->GetPos() << ", " << my_shaftA->GetPosDt() << ", " + << my_shaftA->GetPosDt2() << ", " << my_shaftB->GetPos() << ", " << my_shaftB->GetPosDt() << ", " + << my_shaftB->GetPosDt2() << ", " << my_shaft_gearAB->GetReaction1() << ", " + << my_shaft_gearAB->GetReaction2() << std::endl; } - if (false) { - // - // EXAMPLE 1: - // - - std::cout << "\n\nExample: create a simple power train with ChShaft objects:" << std::endl; - - // We will model a very basic powertrain with two shafts A and B, - // connected by a reducer [ t ] with transmission ratio 't'. Shafts are - // free to rotate, shaft A has an applied torque Ta, so A and B will - // constantly accelerate. Each shaft must have some inertia, it's like a - // flywheel, marked as || in the following scheme: - // - // A B - // Ta ||---[ t ]---|| - // - - // The physical system: it contains all physical objects. - ChSystemNSC sys; - - // Create a 1-degree-of-freedom '1D' mechanical object, that - // is a ChShaft (an item that can oly rotate, with one inertia value - // and maybe one applied torque). The ChShaft objects do not have - // any meaning in 3d: they are just 'building blocks' for making - // power trains as in imput-output black box schemes. - auto my_shaftA = chrono_types::make_shared(); - my_shaftA->SetInertia(10); - my_shaftA->SetAppliedLoad(6); - sys.AddShaft(my_shaftA); - - // Create another shaft. Note that we use shared pointers for ChShaft - // objects, as we did for ChBody objects. Also, note that we must add them - // to the ChSystemNSC. - auto my_shaftB = chrono_types::make_shared(); - my_shaftB->SetInertia(100); - my_shaftB->SetFixed(false); - sys.AddShaft(my_shaftB); - - // Create a ChShaftsGear, that represents a simplified model - // of a reducer, with transmission ratio t, between two ChShaft objects. - // (Note that you could also build a 3D powertrain by creating full rigid bodies - // of ChBody type and join them using ChLinkLockRevolute, ChLinkLockGear 3D constraints, - // but this would introduce many unnecessary degrees of freedom/constraints - // whereas the 1D items of ChShaft type, in this example, make things much simplier). - auto my_shaft_gearAB = chrono_types::make_shared(); - my_shaft_gearAB->Initialize(my_shaftA, my_shaftB); - my_shaft_gearAB->SetTransmissionRatio(-0.1); // ex., a couple of spur gears with 20 and 200 teeth - sys.Add(my_shaft_gearAB); - - std::cout << "\nHere's the system hierarchy:\n" << std::endl; - sys.ShowHierarchy(std::cout); - - // Simulation loop - double end_time = 2.5; - double step_size = 0.01; - - for (double frame_time = 0.05; frame_time < end_time; frame_time += 0.05) { - // Perform simulation up to frame_time - sys.DoFrameDynamics(frame_time, step_size); - - // Print results - std::cout << "Time: " << frame_time << std::endl - << " shaft A rot: " << my_shaftA->GetPos() << " speed: " << my_shaftA->GetPosDt() - << " accel: " << my_shaftA->GetPosDt2() << std::endl - << " shaft B rot: " << my_shaftB->GetPos() << " speed: " << my_shaftB->GetPosDt() - << " accel: " << my_shaftB->GetPosDt2() << std::endl - << " AB gear, torque on A side: " << my_shaft_gearAB->GetReaction1() - << " AB gear, torque on B side: " << my_shaft_gearAB->GetReaction2() << std::endl; + file_results.close(); +} + +void Example2(const std::string& out_dir, ChSolver::Type solver_type, ChTimestepper::Type integrator_type) { + std::cout << "\n-------------------------------------------------------------------------\n"; + std::cout << "Example: a clutch between two shafts\n" << std::endl; + + // We will model a very basic powertrain with two shafts A and B, + // connected by a clutch [ c ]. Shafts (see flywheels || in scheme) + // starts with nonzero speed, and are free to rotate independently + // until the clutch is activated: since activation, they will decelerate + // until they have the same speed. + // + // A B + // Ta ||---[ c ]---|| + // + + // The physical system: it contains all physical objects. + ChSystemNSC sys; + SetChronoSolver(sys, solver_type, integrator_type); + + // Create a ChShaft that starts with nonzero angular velocity + auto my_shaftA = chrono_types::make_shared(); + my_shaftA->SetInertia(0.5); + my_shaftA->SetPosDt(30); + sys.AddShaft(my_shaftA); + + // Create another ChShaft, with opposite initial angular velocity + auto my_shaftB = chrono_types::make_shared(); + my_shaftB->SetInertia(0.6); + my_shaftB->SetPosDt(-10); + sys.AddShaft(my_shaftB); + + // Create a ChShaftsClutch, that represents a simplified model + // of a clutch between two ChShaft objects (something that limits + // the max transmitted torque, up to slippage). + auto my_shaft_clutchAB = chrono_types::make_shared(); + my_shaft_clutchAB->Initialize(my_shaftA, my_shaftB); + my_shaft_clutchAB->SetTorqueLimit(60); + sys.Add(my_shaft_clutchAB); + + // Let's begin the simulation with the clutch disengaged: + my_shaft_clutchAB->SetModulation(0); + + std::cout << "\nSystem hierarchy" << std::endl; + sys.ShowHierarchy(std::cout); + + std::ofstream file_results(out_dir + "/example2.txt"); + + // Simulation loop + double end_time = 2.5; + double step_size = 0.01; + + for (double frame_time = 0.05; frame_time < end_time; frame_time += 0.05) { + // Perform simulation up to frame_time + sys.DoFrameDynamics(frame_time, step_size); + + // Activate the clutch only after 0.8 seconds of simulation + if (frame_time > 0.8) { + my_shaft_clutchAB->SetModulation(1); } + + // Print results + std::cout << "Time: " << frame_time << std::endl + << " shaft A rot: " << my_shaftA->GetPos() << " speed: " << my_shaftA->GetPosDt() + << " accel: " << my_shaftA->GetPosDt2() << std::endl + << " shaft B rot: " << my_shaftB->GetPos() << " speed: " << my_shaftB->GetPosDt() + << " accel: " << my_shaftB->GetPosDt2() << std::endl + << " AB clutch, torque on A side: " << my_shaft_clutchAB->GetReaction1() + << " AB clutch, torque on B side: " << my_shaft_clutchAB->GetReaction2() << std::endl; + file_results << frame_time << ", " << my_shaftA->GetPos() << ", " << my_shaftA->GetPosDt() << ", " + << my_shaftA->GetPosDt2() << ", " << my_shaftB->GetPos() << ", " << my_shaftB->GetPosDt() << ", " + << my_shaftB->GetPosDt2() << ", " << my_shaft_clutchAB->GetReaction1() << ", " + << my_shaft_clutchAB->GetReaction2() << std::endl; } - if (false) { - // - // EXAMPLE 2: - // - - std::cout << "\n\nExample: a clutch between two shafts" << std::endl; - - // We will model a very basic powertrain with two shafts A and B, - // connected by a clutch [ c ]. Shafts (see flywheels || in scheme) - // starts with nonzero speed, and are free to rotate independently - // until the clutch is activated: since activation, they will decelerate - // until they have the same speed. - // - // A B - // Ta ||---[ c ]---|| - // - - // The physical system: it contains all physical objects. - ChSystemNSC sys; - - // Create a ChShaft that starts with nonzero angular velocity - auto my_shaftA = chrono_types::make_shared(); - my_shaftA->SetInertia(0.5); - my_shaftA->SetPosDt(30); - sys.AddShaft(my_shaftA); - - // Create another ChShaft, with opposite initial angular velocity - auto my_shaftB = chrono_types::make_shared(); - my_shaftB->SetInertia(0.6); - my_shaftB->SetPosDt(-10); - sys.AddShaft(my_shaftB); - - // Create a ChShaftsClutch, that represents a simplified model - // of a clutch between two ChShaft objects (something that limits - // the max transmitted torque, up to slippage). - auto my_shaft_clutchAB = chrono_types::make_shared(); - my_shaft_clutchAB->Initialize(my_shaftA, my_shaftB); - my_shaft_clutchAB->SetTorqueLimit(60); - sys.Add(my_shaft_clutchAB); - - // Let's begin the simulation with the clutch disengaged: - my_shaft_clutchAB->SetModulation(0); - - std::cout << "\nHere's the system hierarchy:\n" << std::endl; - sys.ShowHierarchy(std::cout); - - // Simulation loop - double end_time = 2.5; - double step_size = 0.01; - - for (double frame_time = 0.05; frame_time < end_time; frame_time += 0.05) { - // Perform simulation up to frame_time - sys.DoFrameDynamics(frame_time, step_size); - - // Activate the clutch only after 0.8 seconds of simulation - if (frame_time > 0.8) { - my_shaft_clutchAB->SetModulation(1); - } - - // Print results - std::cout << "Time: " << frame_time << std::endl - << " shaft A rot: " << my_shaftA->GetPos() << " speed: " << my_shaftA->GetPosDt() - << " accel: " << my_shaftA->GetPosDt2() << std::endl - << " shaft B rot: " << my_shaftB->GetPos() << " speed: " << my_shaftB->GetPosDt() - << " accel: " << my_shaftB->GetPosDt2() << std::endl - << " AB clutch, torque on A side: " << my_shaft_clutchAB->GetReaction1() - << " AB clutch, torque on B side: " << my_shaft_clutchAB->GetReaction2() << std::endl; - } + file_results.close(); +} + +void Example3(const std::string& out_dir, ChSolver::Type solver_type, ChTimestepper::Type integrator_type) { + std::cout << "\n-------------------------------------------------------------------------\n"; + std::cout << "Example: an epicycloidal reducer\n" << std::endl; + + // We will model an epicycloidal reducer using the ChShaftsPlanetary + // constraint. + // The ChShaftsPlanetary makes a kinematic constraint between three + // shafts: so one of them will be 'fixed' and will represent the truss + // of the reducer -in s reducer, this is the role of the + // large gear with inner teeth- and the two remaining shafts are the + // input and output shafts (in other cases, such as the differential + // planetary gear of the cars, all three shafts are free). + // Also, a brake is applied for the output shaft: the ChShaftsClutch + // will be used to this end, it's enough that one of the two shafts is fixed. + // In the following scheme, the brake is [ b ], the planetary (the + // reducer) is [ p ], the shafts are A,B,C,D applied torque is Ta, inertias + // of free shafts are shown as flywheels || , and fixed shafts are marked with * . + // + // A B D + // Ta ||---[ p ]---||---[ b ]---* + // [ ]---* + // C + + // The physical system: it contains all physical objects. + ChSystemNSC sys; + SetChronoSolver(sys, solver_type, integrator_type); + + // Create shaft A, with applied torque + auto my_shaftA = chrono_types::make_shared(); + my_shaftA->SetInertia(0.5); + my_shaftA->SetAppliedLoad(10); + sys.AddShaft(my_shaftA); + + // Create shaft B + auto my_shaftB = chrono_types::make_shared(); + my_shaftB->SetInertia(0.5); + sys.AddShaft(my_shaftB); + + // Create shaft C, that will be fixed (to be used as truss of epicycloidal reducer) + auto my_shaftC = chrono_types::make_shared(); + my_shaftC->SetFixed(true); + sys.AddShaft(my_shaftC); + + // Create a ChShaftsPlanetary, that represents a simplified model + // of a planetary gear between THREE ChShaft objects (ex.: a car differential) + // An epicycloidal reducer is a special type of planetary gear. + auto my_shaft_planetaryBAC = chrono_types::make_shared(); + my_shaft_planetaryBAC->Initialize(my_shaftB, my_shaftA, my_shaftC); // output, carrier, fixed + // We can set the ratios of the planetary using a simplified formula, for the + // so called 'Willis' case. Imagine we hold fixed the carrier (shaft B in epic. reducers), + // and leave free the truss C (the outer gear with inner teeth in our reducer); which is + // the transmission ratio t0 that we get? It is simply t0=-Za/Zc, with Z = num of teeth of gears. + // So just use the following to set all the three ratios automatically: + double t0 = -50.0 / 100.0; // suppose, in the reducer, that pinion A has 50 teeth and truss has 100 inner teeth. + my_shaft_planetaryBAC->SetTransmissionRatioOrdinary(t0); + sys.Add(my_shaft_planetaryBAC); + + // Now, let's make a shaft D, that is fixed, and used for the right side + // of a clutch (so the clutch will act as a brake). + auto my_shaftD = chrono_types::make_shared(); + my_shaftD->SetFixed(true); + sys.Add(my_shaftD); + + // Make the brake. It is, in fact a clutch between shafts B and D, where + // D is fixed as a truss, so the clutch will operate as a brake. + auto my_shaft_clutchBD = chrono_types::make_shared(); + my_shaft_clutchBD->Initialize(my_shaftB, my_shaftD); + my_shaft_clutchBD->SetTorqueLimit(60); + sys.Add(my_shaft_clutchBD); + + std::cout << "\nSystem hierarchy" << std::endl; + sys.ShowHierarchy(std::cout); + + std::ofstream file_results(out_dir + "/example3.txt"); + + // Simulation loop + double end_time = 2.5; + double step_size = 0.01; + + for (double frame_time = 0.05; frame_time < end_time; frame_time += 0.05) { + // Perform simulation up to frame_time + sys.DoFrameDynamics(frame_time, step_size); + + // Print results + std::cout << "Time: " << frame_time << std::endl + << " shaft A rot: " << my_shaftA->GetPos() << " speed: " << my_shaftA->GetPosDt() + << " accel: " << my_shaftA->GetPosDt2() << std::endl + << " shaft B rot: " << my_shaftB->GetPos() << " speed: " << my_shaftB->GetPosDt() + << " accel: " << my_shaftB->GetPosDt2() << std::endl + << " epicycloidal react torques on shafts - on A: " << my_shaft_planetaryBAC->GetReaction2() + << " , on B: " << my_shaft_planetaryBAC->GetReaction1() + << " , on C: " << my_shaft_planetaryBAC->GetTorqueReactionOn3() << std::endl; + file_results << frame_time << ", " << my_shaftA->GetPos() << ", " << my_shaftA->GetPosDt() << ", " + << my_shaftA->GetPosDt2() << ", " << my_shaftB->GetPos() << ", " << my_shaftB->GetPosDt() << ", " + << my_shaftB->GetPosDt2() << ", " << my_shaft_planetaryBAC->GetReaction2() << ", " + << my_shaft_planetaryBAC->GetReaction1() << ", " << my_shaft_planetaryBAC->GetTorqueReactionOn3() + << std::endl; } - if (false) { - // - // EXAMPLE 3: - // - - std::cout << "\n\nExample: an epicycloidal reducer" << std::endl; - - // We will model an epicycloidal reducer using the ChShaftsPlanetary - // constraint. - // The ChShaftsPlanetary makes a kinematic constraint between three - // shafts: so one of them will be 'fixed' and will represent the truss - // of the reducer -in s reducer, this is the role of the - // large gear with inner teeth- and the two remaining shafts are the - // input and output shafts (in other cases, such as the differential - // planetary gear of the cars, all three shafts are free). - // Also, a brake is applied for the output shaft: the ChShaftsClutch - // will be used to this end, it's enough that one of the two shafts is fixed. - // In the following scheme, the brake is [ b ], the planetary (the - // reducer) is [ p ], the shafts are A,B,C,D applied torque is Ta, inertias - // of free shafts are shown as flywheels || , and fixed shafts are marked with * . - // - // A B D - // Ta ||---[ p ]---||---[ b ]---* - // [ ]---* - // C - - // The physical system: it contains all physical objects. - ChSystemNSC sys; - - // Create shaft A, with applied torque - auto my_shaftA = chrono_types::make_shared(); - my_shaftA->SetInertia(0.5); - my_shaftA->SetAppliedLoad(10); - sys.AddShaft(my_shaftA); - - // Create shaft B - auto my_shaftB = chrono_types::make_shared(); - my_shaftB->SetInertia(0.5); - sys.AddShaft(my_shaftB); - - // Create shaft C, that will be fixed (to be used as truss of epicycloidal reducer) - auto my_shaftC = chrono_types::make_shared(); - my_shaftC->SetFixed(true); - sys.AddShaft(my_shaftC); - - // Create a ChShaftsPlanetary, that represents a simplified model - // of a planetary gear between THREE ChShaft objects (ex.: a car differential) - // An epicycloidal reducer is a special type of planetary gear. - auto my_shaft_planetaryBAC = chrono_types::make_shared(); - my_shaft_planetaryBAC->Initialize(my_shaftB, my_shaftA, my_shaftC); // output, carrier, fixed - // We can set the ratios of the planetary using a simplified formula, for the - // so called 'Willis' case. Imagine we hold fixed the carrier (shaft B in epic. reducers), - // and leave free the truss C (the outer gear with inner teeth in our reducer); which is - // the transmission ratio t0 that we get? It is simply t0=-Za/Zc, with Z = num of teeth of gears. - // So just use the following to set all the three ratios automatically: - double t0 = - -50.0 / 100.0; // suppose, in the reducer, that pinion A has 50 teeth and truss has 100 inner teeth. - my_shaft_planetaryBAC->SetTransmissionRatioOrdinary(t0); - sys.Add(my_shaft_planetaryBAC); - - // Now, let's make a shaft D, that is fixed, and used for the right side - // of a clutch (so the clutch will act as a brake). - auto my_shaftD = chrono_types::make_shared(); - my_shaftD->SetFixed(true); - sys.Add(my_shaftD); - - // Make the brake. It is, in fact a clutch between shafts B and D, where - // D is fixed as a truss, so the clutch will operate as a brake. - auto my_shaft_clutchBD = chrono_types::make_shared(); - my_shaft_clutchBD->Initialize(my_shaftB, my_shaftD); - my_shaft_clutchBD->SetTorqueLimit(60); - sys.Add(my_shaft_clutchBD); - - std::cout << "\nHere's the system hierarchy:\n" << std::endl; - sys.ShowHierarchy(std::cout); - - // Simulation loop - double end_time = 2.5; - double step_size = 0.01; - - for (double frame_time = 0.05; frame_time < end_time; frame_time += 0.05) { - // Perform simulation up to frame_time - sys.DoFrameDynamics(frame_time, step_size); - - // Print results - std::cout << "Time: " << frame_time << std::endl - << " shaft A rot: " << my_shaftA->GetPos() << " speed: " << my_shaftA->GetPosDt() - << " accel: " << my_shaftA->GetPosDt2() << std::endl - << " shaft B rot: " << my_shaftB->GetPos() << " speed: " << my_shaftB->GetPosDt() - << " accel: " << my_shaftB->GetPosDt2() << std::endl - << " epicycloidal react torques on shafts - on A: " << my_shaft_planetaryBAC->GetReaction2() - << " , on B: " << my_shaft_planetaryBAC->GetReaction1() - << " , on C: " << my_shaft_planetaryBAC->GetTorqueReactionOn3() << std::endl; - } + file_results.close(); +} + +void Example4a(const std::string& out_dir, ChSolver::Type solver_type, ChTimestepper::Type integrator_type) { + std::cout << "\n-------------------------------------------------------------------------\n"; + std::cout << "Example: constraint between a ChBody and a ChShaft\n" << std::endl; + + // Suppose you want to create a 3D model, for instance a slider-crank, + // built with multiple ChBody objects; moreover you want to create a + // powertrain, for instance a motor, a clutch, etc, for the rotation of + // the crank. How to connect the '1D items' of ChShaft class to the 3D + // items of ChBody class? The solution is to use the ChShaftBodyRotation constraint, + // shown as [ bs ] in the following scheme, where the 3D body is shown as <>. + // In this example we also add a 'torsional spring damper' C, shown as [ t ] + // that connects shafts A and C (C is shown as * because fixed). + // + // B A C + // Ta <>---[ bs ]---||---[ t ]---* + // + + // The physical system: it contains all physical objects. + ChSystemNSC sys; + SetChronoSolver(sys, solver_type, integrator_type); + + // Create 'A', a 1D shaft + auto my_shaftA = chrono_types::make_shared(); + my_shaftA->SetInertia(9); + sys.AddShaft(my_shaftA); + + // Create 'C', a 1D shaft, fixed + auto my_shaftC = chrono_types::make_shared(); + my_shaftC->SetFixed(true); + sys.AddShaft(my_shaftC); + + // Create 'B', a 3D rigid body + auto my_bodyB = chrono_types::make_shared(); + my_bodyB->AccumulateTorque(ChVector3d(0, 0, 3), true); // set some constant torque to body + sys.Add(my_bodyB); + + // Make the torsional spring-damper between shafts A and C. + auto my_shaft_torsionAC = chrono_types::make_shared(); + my_shaft_torsionAC->Initialize(my_shaftA, my_shaftC); + my_shaft_torsionAC->SetTorsionalStiffness(40); + my_shaft_torsionAC->SetTorsionalDamping(0); + sys.Add(my_shaft_torsionAC); + + // Make the shaft 'A' connected to the rotation of the 3D body 'B'. + // We must specify the direction (in body coordinates) along which the + // shaft will affect the body. + auto my_shaftbody_connection = chrono_types::make_shared(); + ChVector3d mshaftdir(VECT_Z); + my_shaftbody_connection->Initialize(my_shaftA, my_bodyB, mshaftdir); + sys.Add(my_shaftbody_connection); + + std::cout << "\nSystem hierarchy" << std::endl; + sys.ShowHierarchy(std::cout); + + std::ofstream file_results(out_dir + "/example4a.txt"); + + // Simulation loop + double end_time = 2.5; + double step_size = 0.01; + + for (double frame_time = 0.05; frame_time < end_time; frame_time += 0.05) { + // Perform simulation up to frame_time + sys.DoFrameDynamics(frame_time, step_size); + + // Print results + std::cout << "Time: " << frame_time << std::endl + << " shaft A rot: " << my_shaftA->GetPos() << " speed: " << my_shaftA->GetPosDt() + << " accel: " << my_shaftA->GetPosDt2() << std::endl + << " Body B angular speed on z: " << my_bodyB->GetAngVelLocal().z() + << " accel on z: " << my_bodyB->GetAngAccLocal().z() << std::endl + << " AC spring, torque on A side: " << my_shaft_torsionAC->GetReaction1() + << " torque on C side: " << my_shaft_torsionAC->GetReaction2() << std::endl + << " shafts/body reaction, on shaft A: " << my_shaftbody_connection->GetTorqueReactionOnShaft() + << " , on body (x y z): " << my_shaftbody_connection->GetTorqueReactionOnBody().x() << " " + << my_shaftbody_connection->GetTorqueReactionOnBody().y() << " " + << my_shaftbody_connection->GetTorqueReactionOnBody().z() << " " << std::endl; + file_results << frame_time << ", " << my_shaftA->GetPos() << ", " << my_shaftA->GetPosDt() << ", " + << my_shaftA->GetPosDt2() << ", " << my_bodyB->GetAngVelLocal().z() << ", " + << my_bodyB->GetAngAccLocal().z() << ", " << my_shaft_torsionAC->GetReaction1() << ", " + << my_shaft_torsionAC->GetReaction2() << ", " + << my_shaftbody_connection->GetTorqueReactionOnShaft() << ", " + << my_shaftbody_connection->GetTorqueReactionOnBody().x() << ", " + << my_shaftbody_connection->GetTorqueReactionOnBody().y() << ", " + << my_shaftbody_connection->GetTorqueReactionOnBody().z() << std::endl; } - if (false) { - // - // EXAMPLE 4: - // - - std::cout << "\n\nExample: constraint between a ChBody and a ChShaft" << std::endl; - - // Suppose you want to create a 3D model, for instance a slider-crank, - // built with multiple ChBody objects; moreover you want to create a - // powertrain, for instance a motor, a clutch, etc, for the rotation of - // the crank. How to connect the '1D items' of ChShaft class to the 3D - // items of ChBody class? The solution is to use the ChShaftBodyRotation constraint, - // shown as [ bs ] in the following scheme, where the 3D body is shown as <>. - // In this example we also add a 'torsional spring damper' C, shown as [ t ] - // that connects shafts A and C (C is shown as * because fixed). - // - // B A C - // Ta <>---[ bs ]---||---[ t ]---* - // - - // The physical system: it contains all physical objects. - ChSystemNSC sys; - - // Create 'A', a 1D shaft - auto my_shaftA = chrono_types::make_shared(); - my_shaftA->SetInertia(9); - sys.AddShaft(my_shaftA); - - // Create 'C', a 1D shaft, fixed - auto my_shaftC = chrono_types::make_shared(); - my_shaftC->SetFixed(true); - sys.AddShaft(my_shaftC); - - // Create 'B', a 3D rigid body - auto my_bodyB = chrono_types::make_shared(); - my_bodyB->AccumulateTorque(ChVector3d(0, 0, 3), true); // set some constant torque to body - sys.Add(my_bodyB); - - // Make the torsional spring-damper between shafts A and C. - auto my_shaft_torsionAC = chrono_types::make_shared(); - my_shaft_torsionAC->Initialize(my_shaftA, my_shaftC); - my_shaft_torsionAC->SetTorsionalStiffness(40); - my_shaft_torsionAC->SetTorsionalDamping(0); - sys.Add(my_shaft_torsionAC); - - // Make the shaft 'A' connected to the rotation of the 3D body 'B'. - // We must specify the direction (in body coordinates) along which the - // shaft will affect the body. - auto my_shaftbody_connection = chrono_types::make_shared(); - ChVector3d mshaftdir(VECT_Z); - my_shaftbody_connection->Initialize(my_shaftA, my_bodyB, mshaftdir); - sys.Add(my_shaftbody_connection); - - std::cout << "\nHere's the system hierarchy:\n" << std::endl; - sys.ShowHierarchy(std::cout); - - // Perform a very simple simulation loop.. - double end_time = 2.5; - double step_size = 0.01; - - for (double frame_time = 0.05; frame_time < end_time; frame_time += 0.05) { - // Perform simulation up to frame_time - sys.DoFrameDynamics(frame_time, step_size); - - // Print results - std::cout << "Time: " << frame_time << std::endl - << " shaft A rot: " << my_shaftA->GetPos() << " speed: " << my_shaftA->GetPosDt() - << " accel: " << my_shaftA->GetPosDt2() << std::endl - << " Body B angular speed on z: " << my_bodyB->GetAngVelLocal().z() - << " accel on z: " << my_bodyB->GetAngAccLocal().z() << std::endl - << " AC spring, torque on A side: " << my_shaft_torsionAC->GetReaction1() - << " torque on C side: " << my_shaft_torsionAC->GetReaction2() << std::endl - << " shafts/body reaction, on shaft A: " << my_shaftbody_connection->GetTorqueReactionOnShaft() - << " , on body (x y z): " << my_shaftbody_connection->GetTorqueReactionOnBody().x() << " " - << my_shaftbody_connection->GetTorqueReactionOnBody().y() << " " - << my_shaftbody_connection->GetTorqueReactionOnBody().z() << " " << std::endl; - } + file_results.close(); +} + + +void Example4b(const std::string& out_dir, ChSolver::Type solver_type, ChTimestepper::Type integrator_type) { + std::cout << "\n-------------------------------------------------------------------------\n"; + std::cout << "Example: constraint between a ChBody and a ChShaft\n" << std::endl; + + // Same as Example4a, but with a rotational spring-damper between bodies A and C and a torque applied to shaft B. + // + // B A C + // Ta ||---[ bs ]---<>---[ t ]---* + // + + // The physical system: it contains all physical objects. + ChSystemNSC sys; + SetChronoSolver(sys, solver_type, integrator_type); + + // Create 'B', a 1D shaft + auto my_shaftB = chrono_types::make_shared(); + my_shaftB->SetInertia(9); + my_shaftB->SetAppliedLoad(3); + sys.AddShaft(my_shaftB); + + // Create 'C', a 3D rigid body, fixed + auto my_bodyC = chrono_types::make_shared(); + my_bodyC->SetFixed(true); + sys.AddBody(my_bodyC); + + // Create 'A', a 3D rigid body + auto my_bodyA = chrono_types::make_shared(); + sys.AddBody(my_bodyA); + + // Make a torsional spring between bodies A and C + auto my_rsdaAC = chrono_types::make_shared(); + my_rsdaAC->Initialize(my_bodyA, my_bodyC, ChFramed()); + my_rsdaAC->SetSpringCoefficient(40); + my_rsdaAC->SetDampingCoefficient(0); + sys.Add(my_rsdaAC); + + // Make the shaft 'A' connected to the rotation of the 3D body 'B'. + // We must specify the direction (in body coordinates) along which the shaft will affect the body. + auto my_shaftbody_connection = chrono_types::make_shared(); + ChVector3d mshaftdir(VECT_Z); + my_shaftbody_connection->Initialize(my_shaftB, my_bodyA, mshaftdir); + sys.Add(my_shaftbody_connection); + + std::cout << "\nSystem hierarchy" << std::endl; + sys.ShowHierarchy(std::cout); + + std::ofstream file_results(out_dir + "/example4b.txt"); + + // Simulation loop + double end_time = 2.5; + double step_size = 0.01; + + for (double frame_time = 0.05; frame_time < end_time; frame_time += 0.05) { + // Perform simulation up to frame_time + sys.DoFrameDynamics(frame_time, step_size); + + // Print results + std::cout << "Time: " << frame_time << std::endl + << " shaft B rot: " << my_shaftB->GetPos() << " speed: " << my_shaftB->GetPosDt() + << " accel: " << my_shaftB->GetPosDt2() << std::endl + << " Body A angular speed on z: " << my_bodyA->GetAngVelLocal().z() + << " accel on z: " << my_bodyA->GetAngAccLocal().z() << std::endl + << " AC spring, torque on A side: " << my_rsdaAC->GetReaction1().torque + << " torque on C side: " << my_rsdaAC->GetReaction2().torque << std::endl + << " shafts/body reaction, on shaft B: " << my_shaftbody_connection->GetTorqueReactionOnShaft() + << " , on body A: " << my_shaftbody_connection->GetTorqueReactionOnBody() << std::endl; } - if (true) { - // - // EXAMPLE 5: - // - - std::cout << "\n\nExample 5: torque converter and thermal engine" << std::endl; - - // In this example we use a torque converter. - // The torque converter is represented by a ChShaftsTorqueConverter - // object, that connects three '1D items' of ChShaft class: - // - the input shaft A, ie. the impeller connected to the engine - // - the output shaft B, i.e. the turbine connected to the gears and wheels - // - the stator C, that does not rotate and transmits reaction to the truss. - // In the following scheme, the torque converter is represented as [ tc ], - // and we also add a thermal engine, shown with [ e ], and a breaking torque Tb - // (C is shown as * because fixed). - // - // D A B - // *---[ e ]---||---[ tc ]---|| Tb - // [ ]---* - // C - // - - // The physical system: it contains all physical objects. - ChSystemNSC sys; - - // Create 'A', a 1D shaft - auto my_shaftA = chrono_types::make_shared(); - my_shaftA->SetInertia(1.5); - sys.AddShaft(my_shaftA); - - // Create 'B', a 1D shaft - auto my_shaftB = chrono_types::make_shared(); - my_shaftB->SetInertia(3.2); - my_shaftB->SetAppliedLoad(-5); // apply const braking torque - sys.AddShaft(my_shaftB); - - // Create 'C', a 1D shaft, fixed - auto my_shaftC = chrono_types::make_shared(); - my_shaftC->SetFixed(true); - sys.AddShaft(my_shaftC); - - // Create 'D', a 1D shaft, fixed - auto my_shaftD = chrono_types::make_shared(); - my_shaftD->SetFixed(true); - sys.AddShaft(my_shaftD); - - // Make the torque converter and connect the shafts: - // A (input),B (output), C(truss stator) - auto my_torqueconverter = chrono_types::make_shared(); - my_torqueconverter->Initialize(my_shaftA, my_shaftB, my_shaftC); - sys.Add(my_torqueconverter); - - auto mK = chrono_types::make_shared(); - mK->AddPoint(0.0, 15); - mK->AddPoint(0.25, 15); - mK->AddPoint(0.50, 15); - mK->AddPoint(0.75, 16); - mK->AddPoint(0.90, 18); - mK->AddPoint(1.00, 35); - my_torqueconverter->SetCurveCapacityFactor(mK); - - auto mT = chrono_types::make_shared(); - mT->AddPoint(0.0, 2.00); - mT->AddPoint(0.25, 1.80); - mT->AddPoint(0.50, 1.50); - mT->AddPoint(0.75, 1.15); - mT->AddPoint(1.00, 1.00); - my_torqueconverter->SetCurveTorqueRatio(mT); - - // Create the thermal engine, acting on shaft A, the input to - // the torque converter. Note that the thermal engine also - // requires another shaft D, that is used to transmit the - // reaction torque back to a truss (the motor block). - - // Option 1: use a ChShaftsAppliedTorque, it just applies a torque - // to my_shaftA (say, the crankshaft) and the negative torque - // to my_shaftD (say, the motor block). - // It is a quick approach. But you should take care of changing - // the torque at each timestep if you want to simulate a torque curve... - /* - auto my_motor = chrono_types::make_shared(); - my_motor->Initialize(my_shaftA, my_shaftD); - my_motor->SetTorque(30); - sys.Add(my_motor); - */ - - // Option 2: a more powerful approach where you can - // define a torque curve and a throttle value, using the - // ChShaftsThermalEngine. - - auto my_motor = chrono_types::make_shared(); - my_motor->Initialize(my_shaftA, my_shaftD); - sys.Add(my_motor); - - auto mTw = chrono_types::make_shared(); - mTw->AddPoint(-5, 30); // [rad/s], [Nm] - mTw->AddPoint(0, 30); - mTw->AddPoint(200, 60); - mTw->AddPoint(400, 40); - mTw->AddPoint(450, 0); - mTw->AddPoint(500, -60); // torque curve must be defined beyond max speed too - engine might be 'pulled' - my_motor->SetTorqueCurve(mTw); - - std::cout << "\nHere's the system hierarchy:\n" << std::endl; - sys.ShowHierarchy(std::cout); - - // Simulation loop - double end_time = 2.5; - double step_size = 0.01; - - for (double frame_time = 0.05; frame_time < end_time; frame_time += 0.05) { - // Perform simulation up to frame_time - sys.DoFrameDynamics(frame_time, step_size); - - // Print results - std::cout << "Time: " << frame_time << std::endl - << " shaft A rot: " << my_shaftA->GetPos() << " speed: " << my_shaftA->GetPosDt() - << " accel: " << my_shaftA->GetPosDt2() << std::endl - << " shaft B rot: " << my_shaftB->GetPos() << " speed: " << my_shaftB->GetPosDt() - << " accel: " << my_shaftB->GetPosDt2() << std::endl - << " T.Convert.:" - << " R=" << my_torqueconverter->GetSpeedRatio() - << " Tin=" << my_torqueconverter->GetTorqueReactionOnInput() - << " Tout=" << my_torqueconverter->GetTorqueReactionOnOutput() - << " Tstator=" << my_torqueconverter->GetTorqueReactionOnStator() << std::endl - << " T.Motor: " - << " T(w)=" << my_motor->GetReaction1() << "[Nm]" - << " w=" << my_motor->GetRelativePosDt() << "[rad/s]" << std::endl; - } + file_results.close(); +} + +void Example5(const std::string& out_dir, ChSolver::Type solver_type, ChTimestepper::Type integrator_type) { + std::cout << "\n-------------------------------------------------------------------------\n"; + std::cout << "Example 5: torque converter and thermal engine\n" << std::endl; + + // In this example we use a torque converter. + // The torque converter is represented by a ChShaftsTorqueConverter + // object, that connects three '1D items' of ChShaft class: + // - the input shaft A, ie. the impeller connected to the engine + // - the output shaft B, i.e. the turbine connected to the gears and wheels + // - the stator C, that does not rotate and transmits reaction to the truss. + // In the following scheme, the torque converter is represented as [ tc ], + // and we also add a thermal engine, shown with [ e ], and a breaking torque Tb + // (C is shown as * because fixed). + // + // D A B + // *---[ e ]---||---[ tc ]---|| Tb + // [ ]---* + // C + // + + // The physical system: it contains all physical objects. + ChSystemNSC sys; + SetChronoSolver(sys, solver_type, integrator_type); + + // Create 'A', a 1D shaft + auto my_shaftA = chrono_types::make_shared(); + my_shaftA->SetInertia(1.5); + sys.AddShaft(my_shaftA); + + // Create 'B', a 1D shaft + auto my_shaftB = chrono_types::make_shared(); + my_shaftB->SetInertia(3.2); + my_shaftB->SetAppliedLoad(-5); // apply const braking torque + sys.AddShaft(my_shaftB); + + // Create 'C', a 1D shaft, fixed + auto my_shaftC = chrono_types::make_shared(); + my_shaftC->SetFixed(true); + sys.AddShaft(my_shaftC); + + // Create 'D', a 1D shaft, fixed + auto my_shaftD = chrono_types::make_shared(); + my_shaftD->SetFixed(true); + sys.AddShaft(my_shaftD); + + // Make the torque converter and connect the shafts: + // A (input),B (output), C(truss stator) + auto my_torqueconverter = chrono_types::make_shared(); + my_torqueconverter->Initialize(my_shaftA, my_shaftB, my_shaftC); + sys.Add(my_torqueconverter); + + auto mK = chrono_types::make_shared(); + mK->AddPoint(0.0, 15); + mK->AddPoint(0.25, 15); + mK->AddPoint(0.50, 15); + mK->AddPoint(0.75, 16); + mK->AddPoint(0.90, 18); + mK->AddPoint(1.00, 35); + my_torqueconverter->SetCurveCapacityFactor(mK); + + auto mT = chrono_types::make_shared(); + mT->AddPoint(0.0, 2.00); + mT->AddPoint(0.25, 1.80); + mT->AddPoint(0.50, 1.50); + mT->AddPoint(0.75, 1.15); + mT->AddPoint(1.00, 1.00); + my_torqueconverter->SetCurveTorqueRatio(mT); + + // Create the thermal engine, acting on shaft A, the input to + // the torque converter. Note that the thermal engine also + // requires another shaft D, that is used to transmit the + // reaction torque back to a truss (the motor block). + + // Option 1: use a ChShaftsAppliedTorque, it just applies a torque + // to my_shaftA (say, the crankshaft) and the negative torque + // to my_shaftD (say, the motor block). + // It is a quick approach. But you should take care of changing + // the torque at each timestep if you want to simulate a torque curve... + /* + auto my_motor = chrono_types::make_shared(); + my_motor->Initialize(my_shaftA, my_shaftD); + my_motor->SetTorque(30); + sys.Add(my_motor); + */ + + // Option 2: a more powerful approach where you can + // define a torque curve and a throttle value, using the + // ChShaftsThermalEngine. + + auto my_motor = chrono_types::make_shared(); + my_motor->Initialize(my_shaftA, my_shaftD); + sys.Add(my_motor); + + auto mTw = chrono_types::make_shared(); + mTw->AddPoint(-5, 30); // [rad/s], [Nm] + mTw->AddPoint(0, 30); + mTw->AddPoint(200, 60); + mTw->AddPoint(400, 40); + mTw->AddPoint(450, 0); + mTw->AddPoint(500, -60); // torque curve must be defined beyond max speed too - engine might be 'pulled' + my_motor->SetTorqueCurve(mTw); + + std::cout << "\nSystem hierarchy" << std::endl; + sys.ShowHierarchy(std::cout); + + std::ofstream file_results(out_dir + "/example5.txt"); + + // Simulation loop + double end_time = 2.5; + double step_size = 0.01; + + for (double frame_time = 0.05; frame_time < end_time; frame_time += 0.05) { + // Perform simulation up to frame_time + sys.DoFrameDynamics(frame_time, step_size); + + // Print results + std::cout << "Time: " << frame_time << std::endl + << " shaft A rot: " << my_shaftA->GetPos() << " speed: " << my_shaftA->GetPosDt() + << " accel: " << my_shaftA->GetPosDt2() << std::endl + << " shaft B rot: " << my_shaftB->GetPos() << " speed: " << my_shaftB->GetPosDt() + << " accel: " << my_shaftB->GetPosDt2() << std::endl + << " T.Convert.:" + << " R=" << my_torqueconverter->GetSpeedRatio() + << " Tin=" << my_torqueconverter->GetTorqueReactionOnInput() + << " Tout=" << my_torqueconverter->GetTorqueReactionOnOutput() + << " Tstator=" << my_torqueconverter->GetTorqueReactionOnStator() << std::endl + << " T.Motor: " + << " T(w)=" << my_motor->GetReaction1() << "[Nm]" + << " w=" << my_motor->GetRelativePosDt() << "[rad/s]" << std::endl; + file_results << frame_time << ", " << my_shaftA->GetPos() << my_shaftA->GetPosDt() << ", " + << my_shaftA->GetPosDt2() << ", " << my_shaftB->GetPos() << my_shaftB->GetPosDt() << ", " + << my_shaftB->GetPosDt2() << ", " << my_torqueconverter->GetSpeedRatio() << ", " + << my_torqueconverter->GetTorqueReactionOnInput() << ", " + << my_torqueconverter->GetTorqueReactionOnOutput() << ", " + << my_torqueconverter->GetTorqueReactionOnStator() << ", " << my_motor->GetReaction1() << ", " + << my_motor->GetRelativePosDt() << std::endl; } - if (true) { - // - // EXAMPLE 6: - // - - std::cout << "\n\nExample 6: a ratcheting freewheel, as a one-directional clutch" << std::endl; - - // In this example we use a ratcheting freewheel, that acts as a one-directional - // clutch (where the locking in reverse direction happens only at discrete - // steps, depending on the n.of ratcheting teeths). - // The example consists of: - // - the fixed shaft A - // - the shaft B, that rotates back and forth - // - the shaft C, that rotates only unidirectionally - // - the fixed shaft D - // In the following scheme - // - the freewheel is represented as [ fw ], - // - we also add a motor, shown with [ m ], to generate sinusoidal rotation of B for testing - // - we also add a clutch, shown with [ cl ], just to keep the shaft C "stopped" when not - // pushed by the unidirectional freewheel, otherwise would proceed in one direction forever. - // (A,D are shown as * because fixed). - // - // A B C D - // *---[ m ]---||---[ fw ]---||---[ cl ]---* - // - - // The physical system: it contains all physical objects. - ChSystemNSC sys; - - // Create 'A', a 1D shaft, fixed - auto my_shaftA = chrono_types::make_shared(); - my_shaftA->SetFixed(true); - sys.AddShaft(my_shaftA); - - // Create 'B', a 1D shaft - auto my_shaftB = chrono_types::make_shared(); - my_shaftB->SetInertia(1.5); - sys.AddShaft(my_shaftB); - - // Create 'C', a 1D shaft - auto my_shaftC = chrono_types::make_shared(); - my_shaftC->SetInertia(3.2); - sys.AddShaft(my_shaftC); - - // Create D', a 1D shaft, fixed - auto my_shaftD = chrono_types::make_shared(); - my_shaftD->SetFixed(true); - sys.AddShaft(my_shaftD); - - // Make the motor imposing a test sinusoidal rotation - auto my_motor = chrono_types::make_shared(); - my_motor->Initialize(my_shaftA, my_shaftB); - sys.Add(my_motor); - auto my_sinefunction = chrono_types::make_shared(0.001 + 0.5 * CH_2PI / 20, 1.2); - my_motor->SetPositionFunction(my_sinefunction); - - // Make the freewheel: - auto my_freewheel = chrono_types::make_shared(); - my_freewheel->Initialize(my_shaftB, my_shaftC); - my_freewheel->SetRatchetingModeTeeth(25); - // my_freewheel->SetJammingMode(); // this is like having infinite teeth, i.e. no backlash - // my_freewheel->SetFreeBackward(); // this is to reverse the unidirectional behavior - sys.Add(my_freewheel); - - // Make the clutch that keeps the shaft C in place: - auto my_clutch = chrono_types::make_shared(); - my_clutch->Initialize(my_shaftC, my_shaftD); - my_clutch->SetTorqueLimit(100); - sys.Add(my_clutch); - - std::cout << "\nHere's the system hierarchy:\n" << std::endl; - sys.ShowHierarchy(std::cout); - - std::ofstream file_results(out_dir + "/test_clutch.txt"); - - // Simulation loop - double step = 0.01; - - while (sys.GetChTime() < 5.5) { - // Advance dynamics by one step - sys.DoStepDynamics(step); - - // Print results - std::cout << "Time: " << sys.GetChTime() << std::endl - << " shaft B rot: " << my_shaftB->GetPos() << " speed: " << my_shaftB->GetPosDt() - << " accel: " << my_shaftB->GetPosDt2() << std::endl - << " shaft C rot: " << my_shaftC->GetPos() << " speed: " << my_shaftC->GetPosDt() - << " accel: " << my_shaftC->GetPosDt2() << std::endl - << " Torque: Tmotor=" << my_motor->GetReaction1() - << " Tfreewheel=" << my_freewheel->GetReaction1() << " Tclutch=" << my_clutch->GetReaction1() - << " ratchet vane=" << my_freewheel->GetCurrentTeethVane() << std::endl; - file_results << sys.GetChTime() << ", " << my_shaftB->GetPos() << ", " << my_shaftC->GetPos() << ", " - << my_shaftC->GetPosDt() << ", " << my_clutch->GetReaction1() << ", " - << my_freewheel->GetCurrentTeethVane() << "\n"; - } + file_results.close(); +} + +void Example6(const std::string& out_dir, ChSolver::Type solver_type, ChTimestepper::Type integrator_type) { + std::cout << "\n-------------------------------------------------------------------------\n"; + std::cout << "Example 6: a ratcheting freewheel, as a one-directional clutch\n" << std::endl; + + // In this example we use a ratcheting freewheel, that acts as a one-directional + // clutch (where the locking in reverse direction happens only at discrete + // steps, depending on the n.of ratcheting teeths). + // The example consists of: + // - the fixed shaft A + // - the shaft B, that rotates back and forth + // - the shaft C, that rotates only unidirectionally + // - the fixed shaft D + // In the following scheme + // - the freewheel is represented as [ fw ], + // - we also add a motor, shown with [ m ], to generate sinusoidal rotation of B for testing + // - we also add a clutch, shown with [ cl ], just to keep the shaft C "stopped" when not + // pushed by the unidirectional freewheel, otherwise would proceed in one direction forever. + // (A,D are shown as * because fixed). + // + // A B C D + // *---[ m ]---||---[ fw ]---||---[ cl ]---* + // + + // The physical system: it contains all physical objects. + ChSystemNSC sys; + SetChronoSolver(sys, solver_type, integrator_type); + + // Create 'A', a 1D shaft, fixed + auto my_shaftA = chrono_types::make_shared(); + my_shaftA->SetFixed(true); + sys.AddShaft(my_shaftA); + + // Create 'B', a 1D shaft + auto my_shaftB = chrono_types::make_shared(); + my_shaftB->SetInertia(1.5); + sys.AddShaft(my_shaftB); + + // Create 'C', a 1D shaft + auto my_shaftC = chrono_types::make_shared(); + my_shaftC->SetInertia(3.2); + sys.AddShaft(my_shaftC); + + // Create D', a 1D shaft, fixed + auto my_shaftD = chrono_types::make_shared(); + my_shaftD->SetFixed(true); + sys.AddShaft(my_shaftD); + + // Make the motor imposing a test sinusoidal rotation + auto my_motor = chrono_types::make_shared(); + my_motor->Initialize(my_shaftA, my_shaftB); + sys.Add(my_motor); + auto my_sinefunction = chrono_types::make_shared(0.001 + 0.5 * CH_2PI / 20, 1.2); + my_motor->SetPositionFunction(my_sinefunction); + + // Make the freewheel: + auto my_freewheel = chrono_types::make_shared(); + my_freewheel->Initialize(my_shaftB, my_shaftC); + my_freewheel->SetRatchetingModeTeeth(25); + // my_freewheel->SetJammingMode(); // this is like having infinite teeth, i.e. no backlash + // my_freewheel->SetFreeBackward(); // this is to reverse the unidirectional behavior + sys.Add(my_freewheel); + + // Make the clutch that keeps the shaft C in place: + auto my_clutch = chrono_types::make_shared(); + my_clutch->Initialize(my_shaftC, my_shaftD); + my_clutch->SetTorqueLimit(100); + sys.Add(my_clutch); + + std::cout << "\nSystem hierarchy" << std::endl; + sys.ShowHierarchy(std::cout); + + std::ofstream file_results(out_dir + "/example6.txt"); + + // Simulation loop + double step = 0.01; + + while (sys.GetChTime() < 5.5) { + // Advance dynamics by one step + sys.DoStepDynamics(step); + + // Print results + std::cout << "Time: " << sys.GetChTime() << std::endl + << " shaft B rot: " << my_shaftB->GetPos() << " speed: " << my_shaftB->GetPosDt() + << " accel: " << my_shaftB->GetPosDt2() << std::endl + << " shaft C rot: " << my_shaftC->GetPos() << " speed: " << my_shaftC->GetPosDt() + << " accel: " << my_shaftC->GetPosDt2() << std::endl + << " Torque: Tmotor=" << my_motor->GetReaction1() << " Tfreewheel=" << my_freewheel->GetReaction1() + << " Tclutch=" << my_clutch->GetReaction1() + << " ratchet vane=" << my_freewheel->GetCurrentTeethVane() << std::endl; + file_results << sys.GetChTime() << ", " << my_shaftB->GetPos() << ", " << my_shaftC->GetPos() << ", " + << my_shaftC->GetPosDt() << ", " << my_clutch->GetReaction1() << ", " + << my_freewheel->GetCurrentTeethVane() << std::endl; } + file_results.close(); +} + +int main(int argc, char* argv[]) { + std::cout << "Copyright (c) 2017 projectchrono.org\nChrono version: " << CHRONO_VERSION << std::endl; + + // Create output directory + std::string out_dir = GetChronoOutputPath() + "DEMO_POWERTRAIN"; + if (!filesystem::create_directory(filesystem::path(out_dir))) { + std::cout << "Error creating directory " << out_dir << std::endl; + return 1; + } + + auto solver_type = ChSolver::Type::BARZILAIBORWEIN; + auto integrator_type = ChTimestepper::Type::EULER_IMPLICIT_LINEARIZED; + + Example1(out_dir, solver_type, integrator_type); + Example2(out_dir, solver_type, integrator_type); + Example3(out_dir, solver_type, integrator_type); + Example4a(out_dir, solver_type, integrator_type); + Example4b(out_dir, solver_type, integrator_type); + Example5(out_dir, solver_type, integrator_type); + Example6(out_dir, solver_type, integrator_type); + return 0; } diff --git a/src/demos/mbs/demo_MBS_motors.cpp b/src/demos/mbs/demo_MBS_motors.cpp index 5304d75ff6..195f7bc01f 100644 --- a/src/demos/mbs/demo_MBS_motors.cpp +++ b/src/demos/mbs/demo_MBS_motors.cpp @@ -19,6 +19,7 @@ #include #include "chrono/physics/ChSystemNSC.h" +#include "chrono/physics/ChSystemSMC.h" #include "chrono/physics/ChBodyEasy.h" #include "chrono/physics/ChLinkMotorRotationAngle.h" #include "chrono/physics/ChLinkMotorRotationSpeed.h" @@ -46,32 +47,39 @@ using namespace chrono::irrlicht; using namespace chrono::vsg3d; #endif +#include "demos/SetChronoSolver.h" + using namespace chrono; +// ============================================================================= + ChVisualSystem::Type vis_type = ChVisualSystem::Type::VSG; +ChContactMethod contact_method = ChContactMethod::NSC; ChCollisionSystem::Type collision_type = ChCollisionSystem::Type::BULLET; +// ============================================================================= + // Shortcut function that creates two bodies (a slider and a guide) in a given position, // just to simplify the creation of multiple linear motors in this demo. void CreateSliderGuide(std::shared_ptr& guide, std::shared_ptr& slider, std::shared_ptr material, - ChSystem& sys, + ChSystem* sys, const ChVector3d mpos) { guide = chrono_types::make_shared(4, 0.3, 0.6, 1000, material); guide->SetPos(mpos); guide->SetFixed(true); - sys.Add(guide); + sys->Add(guide); slider = chrono_types::make_shared(0.4, 0.2, 0.5, 1000, material); slider->SetPos(mpos + ChVector3d(0, 0.3, 0)); slider->GetVisualShape(0)->SetColor(ChColor(0.6f, 0.6f, 0.0f)); - sys.Add(slider); + sys->Add(slider); auto obstacle = chrono_types::make_shared(0.4, 0.4, 0.4, 8000, material); obstacle->SetPos(mpos + ChVector3d(1.5, 0.4, 0)); obstacle->GetVisualShape(0)->SetColor(ChColor(0.2f, 0.2f, 0.2f)); - sys.Add(obstacle); + sys->Add(obstacle); } // Shortcut function that creates two bodies (a stator and a rotor) in a given position, @@ -81,36 +89,45 @@ void CreateSliderGuide(std::shared_ptr& guide, void CreateStatorRotor(std::shared_ptr& stator, std::shared_ptr& rotor, std::shared_ptr material, - ChSystem& sys, + ChSystem* sys, const ChVector3d& mpos) { stator = chrono_types::make_shared(ChAxis::Y, 0.5, 0.1, 1000, material); stator->SetPos(mpos); stator->SetRot(QuatFromAngleX(CH_PI_2)); stator->SetFixed(true); - sys.Add(stator); + sys->Add(stator); rotor = chrono_types::make_shared(1, 0.1, 0.1, 1000, material); rotor->SetPos(mpos + ChVector3d(0.5, 0, -0.15)); rotor->GetVisualShape(0)->SetColor(ChColor(0.6f, 0.6f, 0.0f)); - sys.Add(rotor); + sys->Add(rotor); } int main(int argc, char* argv[]) { std::cout << "Copyright (c) 2017 projectchrono.org\nChrono version: " << CHRONO_VERSION << std::endl; - // Create a ChronoENGINE physical system - ChSystemNSC sys; - sys.SetCollisionSystemType(collision_type); + // Create a Chrono system + ChSystem* sys; + switch (contact_method) { + case ChContactMethod::SMC: + sys = new ChSystemSMC; + break; + case ChContactMethod::NSC: + sys = new ChSystemNSC; + break; + } + sys->SetCollisionSystemType(collision_type); // Contact material shared among all objects - auto material = chrono_types::make_shared(); + ChContactMaterialData material_data; + auto material = material_data.CreateMaterial(contact_method); // Create a floor that is fixed (that is used also to represent the absolute reference) auto floorBody = chrono_types::make_shared(20, 2, 20, 3000, material); floorBody->SetPos(ChVector3d(0, -2, 0)); floorBody->SetFixed(true); floorBody->GetVisualShape(0)->SetTexture(GetChronoDataFile("textures/blue.png")); - sys.Add(floorBody); + sys->Add(floorBody); // In the following we will create different types of motors // - rotational motors: examples A.1, A.2, etc. @@ -144,7 +161,7 @@ int main(int argc, char* argv[]) { stator1, // body B (master) ChFrame<>(positionA1) // motor frame, in abs. coords ); - sys.Add(rotmotor1); + sys->Add(rotmotor1); // Create a ChFunction to be used for the ChLinkMotorRotationSpeed auto mwspeed = @@ -190,7 +207,7 @@ int main(int argc, char* argv[]) { stator2, // body B (master) ChFrame<>(positionA2) // motor frame, in abs. coords ); - sys.Add(rotmotor2); + sys->Add(rotmotor2); // Create a ChFunction to be used for the ChLinkMotorRotationAngle auto msineangle = chrono_types::make_shared(CH_PI, 0.05); @@ -221,7 +238,7 @@ int main(int argc, char* argv[]) { stator3, // body B (master) ChFrame<>(positionA3) // motor frame, in abs. coords ); - sys.Add(rotmotor3); + sys->Add(rotmotor3); // The torque(time) function: auto mtorquetime = chrono_types::make_shared(160, 2); @@ -248,7 +265,7 @@ int main(int argc, char* argv[]) { stator4, // body B (master) ChFrame<>(positionA4) // motor frame, in abs. coords ); - sys.Add(rotmotor4); + sys->Add(rotmotor4); // Implement our custom torque function. // We could use pre-defined ChFunction classes like sine, constant, ramp, etc., @@ -323,7 +340,7 @@ int main(int argc, char* argv[]) { stator5, // body B (master) ChFrame<>(positionA5) // motor frame, in abs. coords ); - sys.Add(rotmotor5); + sys->Add(rotmotor5); // You may want to change the inertia of 'inner' 1d shafts, (each has default 1kg/m^2) // Note: they adds up to 3D inertia when 3D parts rotate about the link shaft. @@ -352,7 +369,7 @@ int main(int argc, char* argv[]) { auto my_shaftA = chrono_types::make_shared(); my_shaftA->SetInertia(0.03); - sys.AddShaft(my_shaftA); + sys->AddShaft(my_shaftA); // Create 'DRIVE', the hi-speed motor model - as a simple example use a 'imposed speed' motor: this // is the equivalent of the ChLinkMotorRotationSpeed, but for 1D elements: @@ -361,7 +378,7 @@ int main(int argc, char* argv[]) { my_drive->Initialize(my_shaftA, // A , the rotor of the drive rotmotor5->GetInnerShaft2() // S2, the stator of the drive ); - sys.Add(my_drive); + sys->Add(my_drive); // Create a speed(time) function, and use it in my_drive: auto my_driveangle = chrono_types::make_shared(25 * CH_2PI); // 25 [rps] = 1500 [rpm] my_drive->SetSpeedFunction(my_driveangle); @@ -377,7 +394,7 @@ int main(int argc, char* argv[]) { rotmotor5->GetInnerShaft1() // S1, the output shaft ); my_reducer->SetTransmissionRatioOrdinary(1.0 / 100.0); // ratio between wR/wA - sys.Add(my_reducer); + sys->Add(my_reducer); // Btw: later, if you want, you can access / plot speeds and // torques for whatever part of the driveline by putting lines like the following @@ -417,7 +434,7 @@ int main(int argc, char* argv[]) { guide1, // body B (master) ChFrame<>(positionB1, Q_ROTATE_Z_TO_X) // motor frame, in abs. coords ); - sys.Add(motor1); + sys->Add(motor1); // Create a ChFunction to be used for the ChLinkMotorLinearPosition auto msine = chrono_types::make_shared(1.6, 0.5); @@ -456,7 +473,7 @@ int main(int argc, char* argv[]) { guide2, // body B (master) ChFrame<>(positionB2, Q_ROTATE_Z_TO_X) // motor frame, in abs. coords ); - sys.Add(motor2); + sys->Add(motor2); // Create a ChFunction to be used for the ChLinkMotorLinearSpeed auto msp = chrono_types::make_shared(1.6 * 0.5 * CH_2PI, 0.5, CH_PI_2); @@ -509,7 +526,7 @@ int main(int argc, char* argv[]) { guide3, // body B (master) ChFrame<>(positionB3, Q_ROTATE_Z_TO_X) // motor frame, in abs. coords ); - sys.Add(motor3); + sys->Add(motor3); // Create a ChFunction to be used for F(t) in ChLinkMotorLinearForce. auto mF = chrono_types::make_shared(200); @@ -541,7 +558,7 @@ int main(int argc, char* argv[]) { guide4, // body B (master) ChFrame<>(positionB4, Q_ROTATE_Z_TO_X) // motor frame, in abs. coords ); - sys.Add(motor4); + sys->Add(motor4); // Create a ChFunction that computes F by a user-defined algorithm, as a callback. // One quick option would be to inherit from the ChFunction base class, and implement the GetVal() @@ -648,7 +665,7 @@ int main(int argc, char* argv[]) { guide5, // body B (master) ChFrame<>(positionB5, Q_ROTATE_Z_TO_X) // motor frame, in abs. coords ); - sys.Add(motor5); + sys->Add(motor5); // You may want to change the inertia of 'inner' 1D shafts, ("translating" shafts: each has default 1kg) // Note: they adds up to 3D inertia when 3D parts translate about the link shaft. @@ -665,13 +682,13 @@ int main(int argc, char* argv[]) { auto my_shaftB = chrono_types::make_shared(); my_shaftB->SetInertia(0.33); // [kg/m^2] - sys.AddShaft(my_shaftB); + sys->AddShaft(my_shaftB); auto my_driveli = chrono_types::make_shared(); my_driveli->Initialize(my_shaftB, // B , the rotor of the drive motor5->GetInnerShaft2Rot() // S2rot, the stator of the drive ); - sys.Add(my_driveli); + sys->Add(my_driveli); // Create a angle(time) function. It could be something as simple as // auto my_functangle = chrono_types::make_shared(0, 180); @@ -696,7 +713,7 @@ int main(int argc, char* argv[]) { motor5->GetInnerShaft1Lin() // S1lin, the output shaft ); my_rackpinion->SetTransmissionRatios(-1, -1.0 / 100.0, 1); - sys.Add(my_rackpinion); + sys->Add(my_rackpinion); // Btw: later, if you want, you can access / plot speeds and // torques for whatever part of the driveline by putting lines like the following @@ -744,7 +761,7 @@ int main(int argc, char* argv[]) { guide6, // body B (master) ChFrame<>(positionB6, Q_ROTATE_Z_TO_X) // motor frame, in abs. coords ); - sys.Add(motor6); + sys->Add(motor6); // Create a ChFunction to be used for the ChLinkMotorLinearPosition; // Note! look later in the while{...} simulation loop, we'll continuously @@ -768,7 +785,7 @@ int main(int argc, char* argv[]) { case ChVisualSystem::Type::IRRLICHT: { #ifdef CHRONO_IRRLICHT auto vis_irr = chrono_types::make_shared(); - vis_irr->AttachSystem(&sys); + vis_irr->AttachSystem(sys); vis_irr->SetWindowSize(800, 600); vis_irr->SetWindowTitle("Motors"); vis_irr->Initialize(); @@ -788,7 +805,7 @@ int main(int argc, char* argv[]) { case ChVisualSystem::Type::VSG: { #ifdef CHRONO_VSG auto vis_vsg = chrono_types::make_shared(); - vis_vsg->AttachSystem(&sys); + vis_vsg->AttachSystem(sys); vis_vsg->SetWindowTitle("Motors"); vis_vsg->AddCamera(ChVector3d(4.5, 4.5, -10.5)); vis_vsg->SetWindowSize(ChVector2i(800, 600)); @@ -808,11 +825,12 @@ int main(int argc, char* argv[]) { } } - // Modify some setting of the physical system for the simulation - sys.SetSolverType(ChSolver::Type::PSOR); - sys.GetSolver()->AsIterative()->SetMaxIterations(50); + // Set solver and integrator + auto solver_type = ChSolver::Type::BARZILAIBORWEIN; + auto integrator_type = ChTimestepper::Type::EULER_IMPLICIT_LINEARIZED; + SetChronoSolver(*sys, solver_type, integrator_type); - double timestep = 0.005; + double timestep = 0.001; ChRealtimeStepTimer realtime_timer; while (vis->Run()) { vis->BeginScene(); @@ -821,11 +839,11 @@ int main(int argc, char* argv[]) { // Example B.6 requires the setpoint to be changed in the simulation loop. // For example, use a clamped sinusoidal - double t = sys.GetChTime(); + double t = sys->GetChTime(); double Sp = std::min(std::max(2.6 * std::sin(t * 1.8), -1.4), 1.4); motor6setpoint->SetSetpoint(Sp, t); - sys.DoStepDynamics(timestep); + sys->DoStepDynamics(timestep); realtime_timer.Spin(timestep); }