From 87e5229904aca9743dbfc84a3efb2011624c00fe Mon Sep 17 00:00:00 2001 From: psakievich Date: Tue, 17 Jan 2023 20:37:31 -0700 Subject: [PATCH] F/fsi pitch (#1110) * Add some operators for SixDOF * Add missing copyrights * Rename variables and add unit test * Add hub motion based displacement computations * Minor comments and syntax updates * Add unit test for orientation displacements --- include/aero/aero_utils/displacements.h | 68 +++++++++++++------ include/aero/fsi/FSIUtils.h | 73 +++++++++++++++++++++ include/aero/fsi/FSIturbine.h | 9 +++ src/aero/fsi/FSIturbine.C | 15 ++++- src/aero/fsi/OpenfastFSI.C | 9 +++ unit_tests/aero/CMakeLists.txt | 1 + unit_tests/aero/UnitTestDisplacements.C | 86 +++++++++++++++++++------ unit_tests/aero/UnitTestFSIUtils.C | 71 ++++++++++++++++++++ 8 files changed, 291 insertions(+), 41 deletions(-) create mode 100644 include/aero/fsi/FSIUtils.h create mode 100644 unit_tests/aero/UnitTestFSIUtils.C diff --git a/include/aero/aero_utils/displacements.h b/include/aero/aero_utils/displacements.h index 2ff87fb72..b7ac3b3d2 100644 --- a/include/aero/aero_utils/displacements.h +++ b/include/aero/aero_utils/displacements.h @@ -20,29 +20,49 @@ namespace aero { //! WienerMilenkovic parameter struct SixDOF { - // Kind of dangeraous constructor + // Kind of dangerous constructor SixDOF(double* vec) - : translation_({vec[0], vec[1], vec[2]}), - rotation_({vec[3], vec[4], vec[5]}) + : position_({vec[0], vec[1], vec[2]}), + orientation_({vec[3], vec[4], vec[5]}) { } + SixDOF(vs::Vector transDisp, vs::Vector rotDisp) - : translation_(transDisp), rotation_(rotDisp) + : position_(transDisp), orientation_(rotDisp) { } - vs::Vector translation_; - vs::Vector rotation_; + + vs::Vector position_; + vs::Vector orientation_; }; +KOKKOS_FORCEINLINE_FUNCTION +SixDOF +operator+(const SixDOF& a, const SixDOF& b) +{ + // adding b to a, so pushing b wmp onto a stack + return SixDOF( + a.position_ + b.position_, wmp::push(b.orientation_, a.orientation_)); +} + +KOKKOS_FORCEINLINE_FUNCTION +SixDOF +operator-(const SixDOF& a, const SixDOF& b) +{ + // subtracting b from a, so popping b from the a stack + return SixDOF( + a.position_ - b.position_, wmp::pop(b.orientation_, a.orientation_)); +} + KOKKOS_FORCEINLINE_FUNCTION SixDOF linear_interp_total_displacement( const SixDOF start, const SixDOF end, const double interpFactor) { auto transDisp = wmp::linear_interp_translation( - start.translation_, end.translation_, interpFactor); - auto rotDisp = - wmp::linear_interp_rotation(start.rotation_, end.rotation_, interpFactor); + start.position_, end.position_, interpFactor); + auto rotDisp = wmp::linear_interp_rotation( + start.orientation_, end.orientation_, interpFactor); return SixDOF(transDisp, rotDisp); } @@ -52,9 +72,9 @@ linear_interp_total_velocity( const SixDOF start, const SixDOF end, const double interpFactor) { auto transDisp = wmp::linear_interp_translation( - start.translation_, end.translation_, interpFactor); + start.position_, end.position_, interpFactor); auto rotDisp = wmp::linear_interp_translation( - start.translation_, end.translation_, interpFactor); + start.position_, end.position_, interpFactor); return SixDOF(transDisp, rotDisp); } @@ -65,8 +85,20 @@ vs::Vector local_aero_coordinates( const vs::Vector inertialPos, const SixDOF aeroRefPosition) { - const auto shift = inertialPos - aeroRefPosition.translation_; - return wmp::rotate(aeroRefPosition.rotation_, shift); + const auto shift = inertialPos - aeroRefPosition.position_; + return wmp::rotate(aeroRefPosition.orientation_, shift); +} + +// TODO(psakiev) test this +//! Translate coordinate system for SixDOF variable from inertial to a reference +//! coordinate system +KOKKOS_FORCEINLINE_FUNCTION +SixDOF +local_aero_transformation(const SixDOF& inertialPos, const SixDOF& refPos) +{ + return SixDOF( + local_aero_coordinates(inertialPos.position_, refPos), + wmp::push(refPos.orientation_, inertialPos.orientation_)); } //! Convert one array of 6 deflections (transX, transY, transZ, wmX, wmY, @@ -78,12 +110,12 @@ compute_translational_displacements( const SixDOF deflections, const SixDOF referencePos, const vs::Vector cfdPos) { const auto localPos = local_aero_coordinates(cfdPos, referencePos); - const auto delta = cfdPos - referencePos.translation_; + const auto delta = cfdPos - referencePos.position_; // deflection roations need to be applied from the aerodynamic local frame of // reference const vs::Vector rotation = - wmp::rotate(deflections.rotation_, localPos, true); - return deflections.translation_ + rotation - delta; + wmp::rotate(deflections.orientation_, localPos, true); + return deflections.position_ + rotation - delta; } // TODO(psakiev) this function is a place holder for when we need to add pitch @@ -117,8 +149,8 @@ compute_mesh_velocity( const vs::Vector cfdPos) { const auto pointLocal = local_aero_coordinates(cfdPos, referencePos); - const auto pointRotate = wmp::rotate(totalDis.rotation_, pointLocal); - return totalVel.translation_ + (totalVel.rotation_ ^ pointRotate); + const auto pointRotate = wmp::rotate(totalDis.orientation_, pointLocal); + return totalVel.position_ + (totalVel.orientation_ ^ pointRotate); } } // namespace aero diff --git a/include/aero/fsi/FSIUtils.h b/include/aero/fsi/FSIUtils.h new file mode 100644 index 000000000..49b884892 --- /dev/null +++ b/include/aero/fsi/FSIUtils.h @@ -0,0 +1,73 @@ +// Copyright 2017 National Technology & Engineering Solutions of Sandia, LLC +// (NTESS), National Renewable Energy Laboratory, University of Texas Austin, +// Northwest Research Associates. Under the terms of Contract DE-NA0003525 +// with NTESS, the U.S. Government retains certain rights in this software. +// +// This software is released under the BSD 3-clause license. See LICENSE file +// for more details. +// + +#ifndef FSIUTILS_H +#define FSIUTILS_H +#include "aero/aero_utils/displacements.h" +#include "vs/vector.h" +#include + +namespace fsi { + +//! compute displacements from the net motions of the hub (including rotations) +KOKKOS_FORCEINLINE_FUNCTION +vs::Vector +translation_displacements_from_hub_motion( + const aero::SixDOF& hubRef, + const aero::SixDOF& hubDisp, + const aero::SixDOF& bladeRef) +{ + const auto refRelToHub = + aero::local_aero_coordinates(bladeRef.position_, hubRef); + // apply opposite rotation since it is the displacement i.e. what needs to be + // added to get to the current position + const auto positionDueToRotation = + wmp::rotate(hubDisp.orientation_, refRelToHub, true); + const auto referencePosition = bladeRef.position_ - hubRef.position_; + return positionDueToRotation - referencePosition + hubDisp.position_; +} + +//! assemble the orientation just from changes at the root (no pitch, includes +//! hub) +KOKKOS_FORCEINLINE_FUNCTION vs::Vector +orientation_displacments_from_hub_motion( + const aero::SixDOF& rootRef, + const aero::SixDOF& rootDisp, + const aero::SixDOF& bladeRef) +{ + // subtract reference root orientation from the blade reference orientation, + // and then rotate the blade orientation so it is in the root frame of + // reference + const auto rootRelativeRefOrientation = wmp::rotate( + rootRef.orientation_, + wmp::pop(rootRef.orientation_, bladeRef.orientation_)); + // subtract orientation displacements at the root + const auto rootRelativeDisplacement = + wmp::rotate(rootDisp.orientation_, rootRelativeRefOrientation, true); + // apply displacement in the reference frame first and then rotate out of root + // reference frame + return wmp::push(rootDisp.orientation_, rootRelativeDisplacement); +} + +KOKKOS_FORCEINLINE_FUNCTION +aero::SixDOF +displacements_from_hub_motion( + const aero::SixDOF& hubRef, + const aero::SixDOF& hubDisp, + const aero::SixDOF& rootRef, + const aero::SixDOF& rootDisp, + const aero::SixDOF& bladeRef) +{ + return aero::SixDOF( + translation_displacements_from_hub_motion(hubRef, hubDisp, bladeRef), + orientation_displacments_from_hub_motion(rootRef, rootDisp, bladeRef)); +} + +} // namespace fsi +#endif diff --git a/include/aero/fsi/FSIturbine.h b/include/aero/fsi/FSIturbine.h index 178f4c08c..e83ec7726 100644 --- a/include/aero/fsi/FSIturbine.h +++ b/include/aero/fsi/FSIturbine.h @@ -1,3 +1,12 @@ +// Copyright 2017 National Technology & Engineering Solutions of Sandia, LLC +// (NTESS), National Renewable Energy Laboratory, University of Texas Austin, +// Northwest Research Associates. Under the terms of Contract DE-NA0003525 +// with NTESS, the U.S. Government retains certain rights in this software. +// +// This software is released under the BSD 3-clause license. See LICENSE file +// for more details. +// + #ifndef FSITURBINE_H #define FSITURBINE_H diff --git a/src/aero/fsi/FSIturbine.C b/src/aero/fsi/FSIturbine.C index 65087210e..70f3f218b 100644 --- a/src/aero/fsi/FSIturbine.C +++ b/src/aero/fsi/FSIturbine.C @@ -1,4 +1,14 @@ +// Copyright 2017 National Technology & Engineering Solutions of Sandia, LLC +// (NTESS), National Renewable Energy Laboratory, University of Texas Austin, +// Northwest Research Associates. Under the terms of Contract DE-NA0003525 +// with NTESS, the U.S. Government retains certain rights in this software. +// +// This software is released under the BSD 3-clause license. See LICENSE file +// for more details. +// + #include "aero/fsi/FSIturbine.h" +#include "aero/fsi/FSIUtils.h" #include "utils/ComputeVectorDivergence.h" #include @@ -7,7 +17,6 @@ #include "stk_mesh/base/Field.hpp" #include "master_element/MasterElement.h" #include "master_element/MasterElementFactory.h" -#include "aero/aero_utils/displacements.h" #include "netcdf.h" @@ -1968,8 +1977,8 @@ fsiTurbine::mapDisplacements() } // Now the nacelle - stk::mesh::Selector nacsel(stk::mesh::selectUnion(nacelleParts_)); - const auto& nacbkts = bulk_->get_buckets(stk::topology::NODE_RANK, nacsel); + stk::mesh::Selector nacelle(stk::mesh::selectUnion(nacelleParts_)); + const auto& nacbkts = bulk_->get_buckets(stk::topology::NODE_RANK, nacelle); for (auto b : nacbkts) { for (size_t in = 0; in < b->size(); in++) { auto node = (*b)[in]; diff --git a/src/aero/fsi/OpenfastFSI.C b/src/aero/fsi/OpenfastFSI.C index 522ce114e..2127bd9c1 100644 --- a/src/aero/fsi/OpenfastFSI.C +++ b/src/aero/fsi/OpenfastFSI.C @@ -1,3 +1,12 @@ +// Copyright 2017 National Technology & Engineering Solutions of Sandia, LLC +// (NTESS), National Renewable Energy Laboratory, University of Texas Austin, +// Northwest Research Associates. Under the terms of Contract DE-NA0003525 +// with NTESS, the U.S. Government retains certain rights in this software. +// +// This software is released under the BSD 3-clause license. See LICENSE file +// for more details. +// + #include "aero/fsi/OpenfastFSI.h" #include "aero/fsi/FSIturbine.h" #include diff --git a/unit_tests/aero/CMakeLists.txt b/unit_tests/aero/CMakeLists.txt index a5b416db3..0cf847763 100644 --- a/unit_tests/aero/CMakeLists.txt +++ b/unit_tests/aero/CMakeLists.txt @@ -1,5 +1,6 @@ target_sources(${utest_ex_name} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestWienerMilenkovic.C ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestDisplacements.C + ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestFSIUtils.C ) diff --git a/unit_tests/aero/UnitTestDisplacements.C b/unit_tests/aero/UnitTestDisplacements.C index 30ab4fae3..4d44b92b6 100644 --- a/unit_tests/aero/UnitTestDisplacements.C +++ b/unit_tests/aero/UnitTestDisplacements.C @@ -12,13 +12,28 @@ #include namespace test_displacements { +//! Test that two WM Params give the same end location for a point +void +test_wiener_milenkovic( + vs::Vector goldWmp, vs::Vector testWmp, vs::Vector testPoint, double eps) +{ + auto goldPnt = wmp::rotate(goldWmp, testPoint); + auto testPnt = wmp::rotate(testWmp, testPoint); + EXPECT_NEAR(goldPnt.x(), testPnt.x(), eps) + << "Gold WMP: " << goldWmp << " testWmp: " << testWmp; + EXPECT_NEAR(goldPnt.y(), testPnt.y(), eps) + << "Gold WMP: " << goldWmp << " testWmp: " << testWmp; + EXPECT_NEAR(goldPnt.z(), testPnt.z(), eps) + << "Gold WMP: " << goldWmp << " testWmp: " << testWmp; +} + TEST(AeroDisplacements, creation_from_pointer) { std::vector openfastSurrogate(6, 1.0); aero::SixDOF disp(openfastSurrogate.data()); for (int i = 0; i < 3; ++i) { - EXPECT_DOUBLE_EQ(openfastSurrogate[i], disp.translation_[i]); - EXPECT_DOUBLE_EQ(openfastSurrogate[i + 3], disp.rotation_[i]); + EXPECT_DOUBLE_EQ(openfastSurrogate[i], disp.position_[i]); + EXPECT_DOUBLE_EQ(openfastSurrogate[i + 3], disp.orientation_[i]); } } @@ -26,24 +41,55 @@ TEST(AeroDisplacements, creation_from_vs_vector) { aero::SixDOF disp(vs::Vector::one(), 2.0 * vs::Vector::one()); for (int i = 0; i < 3; ++i) { - EXPECT_DOUBLE_EQ(1.0, disp.translation_[i]); - EXPECT_DOUBLE_EQ(2.0, disp.rotation_[i]); + EXPECT_DOUBLE_EQ(1.0, disp.position_[i]); + EXPECT_DOUBLE_EQ(2.0, disp.orientation_[i]); } } -//! Test that two WM Params give the same end location for a point -void -test_wiener_milenkovic( - vs::Vector goldWmp, vs::Vector testWmp, vs::Vector testPoint, double eps) +TEST(AeroDisplacements, add_six_dof_together) { - auto goldPnt = wmp::rotate(goldWmp, testPoint); - auto testPnt = wmp::rotate(testWmp, testPoint); - EXPECT_NEAR(goldPnt.x(), testPnt.x(), eps) - << "Gold WMP: " << goldWmp << " testWmp: " << testWmp; - EXPECT_NEAR(goldPnt.y(), testPnt.y(), eps) - << "Gold WMP: " << goldWmp << " testWmp: " << testWmp; - EXPECT_NEAR(goldPnt.z(), testPnt.z(), eps) - << "Gold WMP: " << goldWmp << " testWmp: " << testWmp; + const auto dispA = vs::Vector::ihat(); + const auto dispB = vs::Vector::khat(); + + const auto orientA = + wmp::create_wm_param(vs::Vector::jhat(), utils::radians(15.0)); + const auto orientB = + wmp::create_wm_param(vs::Vector::ihat(), utils::radians(15.0)); + + const aero::SixDOF a(dispA, orientA); + const aero::SixDOF b(dispB, orientB); + + const auto c = a + b; + + const auto goldTrans = dispA + dispB; + // adding b to a, so pushing b wmp onto a stack + const auto goldRot = wmp::push(orientB, orientA); + for (int i = 0; i < 3; ++i) { + EXPECT_DOUBLE_EQ(goldTrans[i], c.position_[i]) << i; + } + test_wiener_milenkovic(goldRot, c.orientation_, vs::Vector::one(), 1e-12); +} + +TEST(AeroDisplacements, subtract_six_dof) +{ + const auto dispA = vs::Vector::ihat(); + const auto dispB = vs::Vector::khat(); + + const auto orientA = + wmp::create_wm_param(vs::Vector::jhat(), utils::radians(15.0)); + const auto orientB = + wmp::create_wm_param(vs::Vector::ihat(), utils::radians(15.0)); + + const aero::SixDOF a(dispA, orientA); + const aero::SixDOF b(dispB, orientB); + + const auto c = a - b; + const auto goldTrans = dispA - dispB; + const auto goldRot = wmp::pop(orientB, orientA); + for (int i = 0; i < 3; ++i) { + EXPECT_DOUBLE_EQ(goldTrans[i], c.position_[i]) << i; + } + test_wiener_milenkovic(goldRot, c.orientation_, vs::Vector::one(), 1e-12); } // Test displacements interpolation along the (1,0,0) axis, and for rotations @@ -62,7 +108,7 @@ TEST(AeroDisplacements, linear_interp_total_displacements) aero::linear_interp_total_displacement(start, end, interpFactor); for (int i = 0; i < 3; ++i) { - EXPECT_DOUBLE_EQ(interpFactor, interpDisp.translation_[i]) + EXPECT_DOUBLE_EQ(interpFactor, interpDisp.position_[i]) << "Failed i: " << i; } @@ -71,7 +117,7 @@ TEST(AeroDisplacements, linear_interp_total_displacements) const auto axisOffset = (vs::Vector::one() - axis) * 1e-3; const vs::Vector testPoint = axis * interpFactor + axisOffset; auto wmpGold = wmp::create_wm_param(axis, goldAngle); - test_wiener_milenkovic(wmpGold, interpDisp.rotation_, testPoint, 1e-10); + test_wiener_milenkovic(wmpGold, interpDisp.orientation_, testPoint, 1e-10); } TEST( @@ -92,7 +138,7 @@ TEST( // CFD Pos is using "coordinates" field so it is likely fixed in time. Need to // confirm with Ganesh for now we will treat this as an offset from the // referencePos in openfast - const vs::Vector cfdPos = referencePos.translation_ + vs::Vector::ihat(); + const vs::Vector cfdPos = referencePos.position_ + vs::Vector::ihat(); const aero::SixDOF deflections(vs::Vector::one() * delta, vs::Vector::zero()); @@ -154,7 +200,7 @@ TEST( // CFD Pos is using "coordinates" field so it is likely fixed in time. Need to // confirm with Ganesh for now we will treat this as an offset from the // referencePos in openfast - const vs::Vector cfdPos = referencePos.translation_ + vs::Vector::ihat(); + const vs::Vector cfdPos = referencePos.position_ + vs::Vector::ihat(); const aero::SixDOF deflections( vs::Vector::zero(), wmp::create_wm_param(vs::Vector::jhat(), angleDef)); diff --git a/unit_tests/aero/UnitTestFSIUtils.C b/unit_tests/aero/UnitTestFSIUtils.C new file mode 100644 index 000000000..a81086adc --- /dev/null +++ b/unit_tests/aero/UnitTestFSIUtils.C @@ -0,0 +1,71 @@ +// Copyright 2017 National Technology & Engineering Solutions of Sandia, LLC +// (NTESS), National Renewable Energy Laboratory, University of Texas Austin, +// Northwest Research Associates. Under the terms of Contract DE-NA0003525 +// with NTESS, the U.S. Government retains certain rights in this software. +// +// This software is released under the BSD 3-clause license. See LICENSE file +// for more details. + +#include "gtest/gtest.h" +#include "aero/fsi/FSIUtils.h" + +namespace aero { +namespace { + +TEST(FSIUtils, translation_component_from_blade_rotation_displacement) +{ + // hub at (1,1,1), orientation is such that upstream is -ihat + // (upstream reference hub for this test is considered local ihat) + const SixDOF hubReference( + vs::Vector::one(), wmp::create_wm_param(vs::Vector::khat(), M_PI)); + + // blade point of reference is (1,1,2), and twisted 45 degrees relative to + // hub + const SixDOF bladeReference( + vs::Vector(1.0, 1.0, 2.0), + wmp::create_wm_param(vs::Vector::khat(), M_PI_4)); + + // hub has rotated 90 degrees about the upstream (-ihat) direction + const SixDOF hubDisp( + vs::Vector::zero(), wmp::create_wm_param(vs::Vector::ihat(-1.0), M_PI_2)); + + const vs::Vector newPosition(1.0, 0.0, 1.0); + + const auto goldRotationDisp = newPosition - bladeReference.position_; + + const auto bldCurRelativeToHub = + fsi::translation_displacements_from_hub_motion( + hubReference, hubDisp, bladeReference); + + EXPECT_DOUBLE_EQ(goldRotationDisp.x(), bldCurRelativeToHub.x()); + EXPECT_DOUBLE_EQ(goldRotationDisp.y(), bldCurRelativeToHub.y()); + EXPECT_DOUBLE_EQ(goldRotationDisp.z(), bldCurRelativeToHub.z()); +} + +TEST(FSIUtils, orientation_component_from_blade_rotation_displacement) +{ + const auto rootRefOrient = wmp::create_wm_param(vs::Vector::khat(), M_PI); + // blade point oriented 90 degrees off from the root + const auto bladeRefOrient = wmp::create_wm_param(vs::Vector::khat(), M_PI_2); + // net rotation of 90 degrees about the hub + const auto rootDispOrient = wmp::create_wm_param(vs::Vector::ihat(), M_PI_2); + + const SixDOF rootRef(vs::Vector::zero(), rootRefOrient); + const SixDOF bladeRef(vs::Vector::zero(), bladeRefOrient); + const SixDOF rootDisp(vs::Vector::zero(), rootDispOrient); + + const auto netOrientation = + fsi::orientation_displacments_from_hub_motion(rootRef, rootDisp, bladeRef); + + const auto testPoint = vs::Vector::one(); + // rotate 90 about khat then 90 about ihat + const vs::Vector goldEndLocation(-1.0, -1.0, 1.0); + + const auto endLocation = wmp::rotate(netOrientation, testPoint); + + EXPECT_DOUBLE_EQ(goldEndLocation.x(), endLocation.x()); + EXPECT_DOUBLE_EQ(goldEndLocation.y(), endLocation.y()); + EXPECT_DOUBLE_EQ(goldEndLocation.z(), endLocation.z()); +} +} // namespace +} // namespace aero