|
| 1 | +#include <chrono> |
| 2 | +#include <filesystem> |
| 3 | +#include <iomanip> |
| 4 | +#include <vector> |
| 5 | + |
| 6 | +#include <chrono/physics/ChBodyEasy.h> |
| 7 | +#include <chrono/physics/ChSystemNSC.h> |
| 8 | + |
| 9 | +#include <hydroc/helper.h> |
| 10 | +#include <hydroc/hydro_system.h> |
| 11 | + |
| 12 | +using namespace chrono; |
| 13 | + |
| 14 | +int main(int argc, char* argv[]) { |
| 15 | + std::cout << "Chrono version: " << CHRONO_VERSION << "\n\n"; |
| 16 | + |
| 17 | + std::string data_dir; |
| 18 | + if (!hydroc::SetInitialEnvironment(data_dir)) return 1; |
| 19 | + |
| 20 | + std::filesystem::path DATADIR(hydroc::getDataDir()); |
| 21 | + |
| 22 | + auto body1_meshfame = (DATADIR / "demos" / "oswec" / "geometry" / "flap.obj").lexically_normal().generic_string(); |
| 23 | + auto body2_meshfame = (DATADIR / "demos" / "oswec" / "geometry" / "base.obj").lexically_normal().generic_string(); |
| 24 | + auto h5fname = (DATADIR / "demos" / "oswec" / "hydroData" / "oswec.h5").lexically_normal().generic_string(); |
| 25 | + |
| 26 | + ChSystemNSC system; |
| 27 | + system.SetGravitationalAcceleration(ChVector3d(0.0, 0.0, -9.81)); |
| 28 | + double timestep = 0.015; |
| 29 | + system.SetSolverType(ChSolver::Type::GMRES); |
| 30 | + double simulationDuration = hydroc::getSimDuration(200.0, 400.0); |
| 31 | + |
| 32 | + std::vector<double> time_vector; |
| 33 | + std::vector<double> flap_rot; |
| 34 | + |
| 35 | + std::cout << "Attempting to open mesh file: " << body1_meshfame << std::endl; |
| 36 | + std::shared_ptr<ChBody> flap_body = chrono_types::make_shared<ChBodyEasyMesh>( |
| 37 | + body1_meshfame, |
| 38 | + 1000, |
| 39 | + false, |
| 40 | + true, |
| 41 | + false |
| 42 | + ); |
| 43 | + |
| 44 | + system.Add(flap_body); |
| 45 | + flap_body->SetName("body1"); |
| 46 | + flap_body->SetPos(ChVector3d(0, 0, -3.9)); |
| 47 | + flap_body->SetMass(127000.0); |
| 48 | + flap_body->SetInertiaXX(ChVector3d(1.85e6, 1.85e6, 1.85e6)); |
| 49 | + |
| 50 | + std::cout << "Attempting to open mesh file: " << body2_meshfame << std::endl; |
| 51 | + std::shared_ptr<ChBody> base_body = chrono_types::make_shared<ChBodyEasyMesh>( |
| 52 | + body2_meshfame, |
| 53 | + 1000, |
| 54 | + false, |
| 55 | + true, |
| 56 | + false |
| 57 | + ); |
| 58 | + |
| 59 | + system.Add(base_body); |
| 60 | + base_body->SetName("body2"); |
| 61 | + base_body->SetPos(ChVector3d(0, 0, -10.15)); |
| 62 | + base_body->SetMass(999); |
| 63 | + base_body->SetInertiaXX(ChVector3d(1, 1, 1)); |
| 64 | + |
| 65 | + auto ground = chrono_types::make_shared<ChBody>(); |
| 66 | + system.AddBody(ground); |
| 67 | + ground->SetPos(ChVector3d(0, 0, -10.15)); |
| 68 | + ground->SetTag(-1); |
| 69 | + ground->SetFixed(true); |
| 70 | + ground->EnableCollision(false); |
| 71 | + |
| 72 | + auto anchor = chrono_types::make_shared<ChLinkMateGeneric>(); |
| 73 | + anchor->Initialize(base_body, ground, false, base_body->GetVisualModelFrame(), base_body->GetVisualModelFrame()); |
| 74 | + system.Add(anchor); |
| 75 | + anchor->SetConstrainedCoords(true, true, true, true, true, true); |
| 76 | + |
| 77 | + ChQuaternion<> revoluteRot = QuatFromAngleX(CH_PI / 2.0); |
| 78 | + auto revolute = chrono_types::make_shared<ChLinkLockRevolute>(); |
| 79 | + revolute->Initialize(base_body, flap_body, ChFramed(ChVector3d(0.0, 0.0, -8.9), revoluteRot)); |
| 80 | + system.AddLink(revolute); |
| 81 | + |
| 82 | + // PTO damper between flap and base |
| 83 | + auto pto = chrono_types::make_shared<ChLinkRSDA>(); |
| 84 | + pto->Initialize(flap_body, base_body, ChFramed(ChVector3d(0.0, 0.0, -8.9), revoluteRot)); |
| 85 | + pto->SetSpringCoefficient(0.0); |
| 86 | + pto->SetDampingCoefficient(12.0e6); |
| 87 | + system.AddLink(pto); |
| 88 | + |
| 89 | + std::vector<std::shared_ptr<ChBody>> bodies; |
| 90 | + bodies.push_back(flap_body); |
| 91 | + bodies.push_back(base_body); |
| 92 | + |
| 93 | + IrregularWaveParams wave_inputs; |
| 94 | + wave_inputs.ramp_duration = 30.0; |
| 95 | + wave_inputs.wave_height = 2.5; |
| 96 | + wave_inputs.wave_period = 8.0; |
| 97 | + wave_inputs.frequency_min = 0.001; |
| 98 | + wave_inputs.frequency_max = 1.0; |
| 99 | + wave_inputs.nfrequencies = 1000; |
| 100 | + |
| 101 | + std::shared_ptr<IrregularWaves> my_hydro_inputs; |
| 102 | + |
| 103 | + try { |
| 104 | + my_hydro_inputs = std::make_shared<IrregularWaves>(wave_inputs); |
| 105 | + } catch (const std::exception& e) { |
| 106 | + std::cerr << "Caught exception: " << e.what() << '\n'; |
| 107 | + return 1; |
| 108 | + } catch (...) { |
| 109 | + std::cerr << "Caught unknown exception.\n"; |
| 110 | + return 1; |
| 111 | + } |
| 112 | + |
| 113 | + if (!my_hydro_inputs) { |
| 114 | + std::cerr << "ERROR: Failed to create IrregularWaves object." << std::endl; |
| 115 | + return 1; |
| 116 | + } |
| 117 | + |
| 118 | + HydroForces hydro_forces(bodies, h5fname); |
| 119 | + hydro_forces.AddWaves(my_hydro_inputs); |
| 120 | + |
| 121 | + auto start = std::chrono::high_resolution_clock::now(); |
| 122 | + |
| 123 | + while (system.GetChTime() <= simulationDuration) { |
| 124 | + system.DoStepDynamics(timestep); |
| 125 | + time_vector.push_back(system.GetChTime()); |
| 126 | + flap_rot.push_back(flap_body->GetRot().GetCardanAnglesXYZ().y()); |
| 127 | + } |
| 128 | + |
| 129 | + auto end = std::chrono::high_resolution_clock::now(); |
| 130 | + unsigned duration = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); |
| 131 | + std::cout << "Simulation completed in " << duration << " ms" << std::endl; |
| 132 | + |
| 133 | + std::string out_dir = hydroc::getTestOutDir() + "/" + RESULTS_DIR_NAME; |
| 134 | + std::filesystem::create_directories(std::filesystem::path(out_dir)); |
| 135 | + |
| 136 | + std::ofstream outputFile(out_dir + "/" + RESULTS_FILE_NAME + ".txt"); |
| 137 | + if (outputFile.is_open()) { |
| 138 | + outputFile << std::left << std::setw(10) << "Time (s)" << std::right << std::setw(16) |
| 139 | + << "Flap Rotation y (radians)" << std::right << std::setw(16) << "Flap Rotation y (degrees)" |
| 140 | + << std::endl; |
| 141 | + for (size_t i = 0; i < time_vector.size(); ++i) |
| 142 | + outputFile << std::left << std::setw(10) << std::setprecision(3) << std::fixed << time_vector[i] |
| 143 | + << std::right << std::setw(16) << std::setprecision(6) << std::fixed << flap_rot[i] |
| 144 | + << std::right << std::setw(16) << std::setprecision(6) << std::fixed |
| 145 | + << flap_rot[i] * 360.0 / (2.0 * CH_PI) << std::endl; |
| 146 | + outputFile.close(); |
| 147 | + } else { |
| 148 | + std::cout << "Error: Could not open output file for writing." << std::endl; |
| 149 | + return 1; |
| 150 | + } |
| 151 | + |
| 152 | + return 0; |
| 153 | +} |
0 commit comments