Add functions for fluid acceleration from waves

Author: tridelat

include/hydroc/wave_types.h

@@ -70,6 +70,8 @@ class WaveBase {
 
     virtual Eigen::Vector3d GetVelocity(const Eigen::Vector3d& position, double time) = 0;
 
+    virtual Eigen::Vector3d GetAcceleration(const Eigen::Vector3d& position, double time) = 0;
+
     /// @brief Mean water level
     double mwl_ = 0.0;
     /// @brief Gravitational acceleration
@@ -102,6 +104,9 @@ class NoWave : public WaveBase {
     Eigen::Vector3d GetVelocity(const Eigen::Vector3d& position, double time) override {
         return Eigen::Vector3d(0.0, 0.0, 0.0);
     }
+    Eigen::Vector3d GetAcceleration(const Eigen::Vector3d& position, double time) override {
+        return Eigen::Vector3d(0.0, 0.0, 0.0);
+    }
 
   private:
     unsigned int num_bodies_;
@@ -175,6 +180,8 @@ class RegularWave : public WaveBase {
 
     Eigen::Vector3d GetVelocity(const Eigen::Vector3d& position, double time) override;
 
+    Eigen::Vector3d GetAcceleration(const Eigen::Vector3d& position, double time) override;
+
   private:
     unsigned int num_bodies_;
     const WaveMode mode_ = WaveMode::regular;
@@ -365,6 +372,8 @@ class IrregularWaves : public WaveBase {
 
     Eigen::Vector3d GetVelocity(const Eigen::Vector3d& position, double time) override;
 
+    Eigen::Vector3d GetAcceleration(const Eigen::Vector3d& position, double time) override;
+
   private:
     IrregularWaveParams params_;
     std::vector<double> spectrum_;

src/wave_types.cpp

@@ -85,6 +85,37 @@ Eigen::Vector3d GetWaterVelocity(const Eigen::Vector3d& position,
     return water_velocity;
 }
 
+Eigen::Vector3d GetWaterAcceleration(const Eigen::Vector3d& position,
+                                     double time,
+                                     double omega,
+                                     double amplitude,
+                                     double phase,
+                                     double wavenumber,
+                                     double water_depth,
+                                     double mwl) {
+    // assuming wave along global X axis position
+    auto x_pos = position.x();
+    // position relative to mean water level
+    auto z_pos = position.z() - mwl;
+
+    // get water velocity
+    auto water_acceleration = Eigen::Vector3d(0.0, 0.0, 0.0);
+    if (2 * M_PI / wavenumber > water_depth || wavenumber * water_depth > 500.0) {
+        // deep water
+        water_acceleration[0] =
+            omega * omega * amplitude * std::exp(wavenumber * z_pos) * sin(wavenumber * x_pos - omega * time + phase);
+        water_acceleration[2] =
+            -omega * omega * amplitude * std::exp(wavenumber * z_pos) * cos(wavenumber * x_pos - omega * time + phase);
+    } else {
+        // shallow water
+        water_acceleration[0] = omega * omega * amplitude * std::cosh(wavenumber * (z_pos + water_depth)) /
+                                std::sinh(wavenumber * water_depth) * sin(wavenumber * x_pos - omega * time + phase);
+        water_acceleration[2] = -omega * omega * amplitude * std::sinh(wavenumber * (z_pos + water_depth)) /
+                                std::sinh(wavenumber * water_depth) * cos(wavenumber * x_pos - omega * time + phase);
+    }
+    return water_acceleration;
+}
+
 Eigen::Vector3d GetWaterVelocityIrregular(const Eigen::Vector3d& position,
                                           double time,
                                           const Eigen::VectorXd& freqs_hz,
@@ -104,6 +135,25 @@ Eigen::Vector3d GetWaterVelocityIrregular(const Eigen::Vector3d& position,
     return water_velocity;
 }
 
+Eigen::Vector3d GetWaterAccelerationIrregular(const Eigen::Vector3d& position,
+                                              double time,
+                                              const Eigen::VectorXd& freqs_hz,
+                                              const Eigen::VectorXd& spectral_densities,
+                                              const Eigen::VectorXd& spectral_widths,
+                                              const Eigen::VectorXd& wave_phases,
+                                              const Eigen::VectorXd& wavenumbers,
+                                              double water_depth,
+                                              double mwl) {
+    auto water_acceleration = Eigen::Vector3d(0.0, 0.0, 0.0);
+    for (size_t i = 0; i < freqs_hz.size(); ++i) {
+        auto amplitude = std::sqrt(2 * spectral_densities[i] * spectral_widths[i]);
+        auto omega     = 2 * M_PI * freqs_hz[i];
+        water_acceleration +=
+            GetWaterAcceleration(position, time, omega, amplitude, wave_phases[i], wavenumbers[i], water_depth, mwl);
+    }
+    return water_acceleration;
+}
+
 double ComputeWaveNumber(double omega,
                          double water_depth,
                          double g,
@@ -182,6 +232,11 @@ Eigen::Vector3d RegularWave::GetVelocity(const Eigen::Vector3d& position, double
                             wavenumber_, water_depth_, mwl_);
 };
 
+Eigen::Vector3d RegularWave::GetAcceleration(const Eigen::Vector3d& position, double time) {
+    return GetWaterAcceleration(position, time, regular_wave_omega_, regular_wave_amplitude_, regular_wave_phase_,
+                                wavenumber_, water_depth_, mwl_);
+};
+
 double RegularWave::GetElevation(const Eigen::Vector3d& position, double time) {
     return GetEta(position, time, regular_wave_omega_, regular_wave_amplitude_, regular_wave_phase_, wavenumber_);
 };
@@ -400,6 +455,22 @@ Eigen::Vector3d IrregularWaves::GetVelocity(const Eigen::Vector3d& position, dou
                                      spectral_widths_, wave_phases_, wavenumbers_, water_depth_, mwl_);
 };
 
+Eigen::Vector3d IrregularWaves::GetAcceleration(const Eigen::Vector3d& position, double time) {
+    // apply wave stretching (if enabled)
+    auto position_stretched = position;
+    if (params_.wave_stretching_) {
+        auto eta = GetEtaIrregular(position, time, spectrum_frequencies_, spectral_densities_, spectral_widths_,
+                                   wave_phases_, wavenumbers_);
+        // position relative to mean water level
+        auto z_pos = position.z() - mwl_;
+        // Wheeler stretching
+        position_stretched[2] = water_depth_ * (z_pos - eta) / (water_depth_ + eta);
+    }
+
+    return GetWaterAccelerationIrregular(position_stretched, time, spectrum_frequencies_, spectral_densities_,
+                                         spectral_widths_, wave_phases_, wavenumbers_, water_depth_, mwl_);
+};
+
 double IrregularWaves::GetElevation(const Eigen::Vector3d& position, double time) {
     return GetEtaIrregular(position, time, spectrum_frequencies_, spectral_densities_, spectral_widths_, wave_phases_,
                            wavenumbers_);