Adding tidal potential to gravity model

doc/modules/changes/20250716_hyunseong96

@@ -0,0 +1,4 @@
+Added: ASPECT now has a new gravity model plugin called 'Radial linear with tidal potnetial' for tidal forces.
+This plugin is useful for modeling long-term interior and surface evolution in moons orbiting a large planet.
+<br>
+(Hyunseong Kim, Antoniette Greta Grima, Wolfgang Bangerth 2025/07/16)

include/aspect/gravity_model/radial.h

@@ -135,6 +135,52 @@ namespace aspect
         double magnitude_at_bottom;
 
     };
+
+
+
+    template <int dim>
+    class RadialLinearWithTidalPotential : public RadialLinear<dim>
+    {
+      public:
+        /**
+         * Return the gravity vector as a function of position.
+         */
+        Tensor<1,dim> gravity_vector (const Point<dim> &position) const override;
+
+        /**
+         * Declare the parameters this class takes through input files.
+         */
+        static
+        void
+        declare_parameters (ParameterHandler &prm);
+
+        /**
+         * Read the parameters this class declares from the parameter file.
+         */
+        void
+        parse_parameters (ParameterHandler &prm) override;
+
+      private:
+        /**
+         * Mass of hosting planet
+         */
+        double M_p;
+
+        /**
+         * Semimajor axis of satellite's orbit
+         */
+        double a_s;
+
+        /**
+         * Rate of non-synchronous rotation in m/year
+         */
+        double b_NSR;
+
+        /**
+         * Radius of modelling satellite
+         */
+        double R1;
+    };
   }
 }
 

source/gravity_model/radial.cc

@@ -25,6 +25,11 @@
 
 #include <deal.II/base/tensor.h>
 
+#include <aspect/geometry_model/spherical_shell.h>
+#include <aspect/geometry_model/sphere.h>
+#include <aspect/geometry_model/chunk.h>
+#include <aspect/geometry_model/ellipsoidal_chunk.h>
+
 namespace aspect
 {
   namespace GravityModel
@@ -178,6 +183,159 @@ namespace aspect
                                "do not have either a spherical or ellipsoidal natural coordinate system."));
 
     }
+
+
+
+// ----------------------------- RadialLinearWithTidalPotential ----------------------
+
+
+    template <int dim>
+    Tensor<1,dim>
+    RadialLinearWithTidalPotential<dim>::gravity_vector (const Point<dim> &p) const
+    {
+      /**
+       * This gravity model is a plugin for a case where the
+       * tidal potential by flattening and non-synchronnous rotation changes gravity with position and time.
+       *
+       * The equation implemented in this heating model is from Tobie et al. (2025) (https://doi.org/10.1007/s11214-025-01136-y),
+       * which is defined as:
+       * g = -magnitude - gradient (-tidal potential)).
+       * potential = (3 G M_p) / (2 a_s^3) * r^2 * (Tstar + T0)
+       * Tstar = 1/6 *(1-3*cos(theta)^2) and T0=1/2sin(theta)^2*cos(2*lambda + 2*b*t)
+       * where G = gravitational constant, M_p = mass of planet, a_s = semimajor axis of satellite's orbit, b = angular rate of nonsynchronous rotation.
+       * r, theta and lambda are radial distance, polar angle and azimuthal angle, respectively.
+       * b [1/s] = b_NSR [m/yr] / (circumference of satellite [m]) / year_to_seconds [s/yr] * 2 * pi
+       *
+       * Notation of this potential equation is converted from spherical coordinates to cartesian coordinates.
+       * Therefore, gradient of potential is (3 G M_p) / (2 a_s^3) * ( 1 / 6 * ( x^2 + y^2 - 2 * z^2) + 1 / 2 * (C1*(x^2 + y^2) - 2 * C2 * x * y)))
+       * where C1 = cos(2*b*t) and C2 = sin(2*b*t)
+       */
+      const Tensor<1,dim> e_x = Point<dim>::unit_vector(0);
+      const Tensor<1,dim> e_y = Point<dim>::unit_vector(1);
+      const Tensor<1,dim> e_z = Point<dim>::unit_vector(2);
+
+      const double t = (this->simulator_is_past_initialization()) ? this->get_time() : 0.0;
+      const double x = p[0];
+      const double y = p[1];
+      const double z = p[2];
+
+      const double C1 = std::cos( 2. * b_NSR / R1 / year_in_seconds * t);
+      const double C2 = std::sin( 2. * b_NSR / R1 / year_in_seconds * t);
+
+      const double dTstar_over_dx = 1. / 6. * ( 2. * x );
+      const double dT0_over_dx    = 1. / 2. * ( 2. * C1 * x - 2. * C2 * y );
+
+      const double dTstar_over_dy = 1. / 6. * ( 2. * y );
+      const double dT0_over_dy    = 1. / 2. * ( -2. * C1 * y - 2. * C2 * x );
+
+      const double dTstar_over_dz = 1. / 6. * ( -4. * z );
+      const double dT0_over_dz    = 0;
+
+      const double G = aspect::constants::big_g;
+      const double T_factor = 3. * G * M_p / 2. / a_s / a_s / a_s;
+
+      const Tensor<1,dim> tidal_gravity = T_factor *
+                                          ( (dTstar_over_dx + dT0_over_dx) * e_x
+                                            + (dTstar_over_dy + dT0_over_dy) * e_y
+                                            + (dTstar_over_dz + dT0_over_dz) * e_z);
+
+      return RadialLinear<dim>::gravity_vector(p) + tidal_gravity;
+    }
+
+
+    template <int dim>
+    void
+    RadialLinearWithTidalPotential<dim>::declare_parameters (ParameterHandler &prm)
+    {
+      RadialLinear<dim>::declare_parameters (prm);
+
+      prm.enter_subsection("Gravity model");
+      {
+        prm.enter_subsection("Radial linear with tidal potnetial");
+        {
+          prm.declare_entry ("Mass of planet", "1.898e27",
+                             Patterns::Double (),
+                             "Mass of satellte's host planet. "
+                             "Default value is for Europa, therefore, mass of Jupiter. "
+                             "Units is $kg$.");
+          prm.declare_entry ("Semimajor axis of satellite", "670900000",
+                             Patterns::Double (),
+                             "Length of semimajor axis of satellite. "
+                             "Default value is for Europa's semimajor axis"
+                             "Units is $m$.");
+          prm.declare_entry ("Rate of nonsynchronous rotation", "1000",
+                             Patterns::Double (),
+                             "Rate of nonsynchronous rotation (NSR). "
+                             "This will be converted to angular rate. "
+                             "Units is $m/year$");
+        }
+        prm.leave_subsection ();
+      }
+      prm.leave_subsection ();
+    }
+
+
+    template <int dim>
+    void
+    RadialLinearWithTidalPotential<dim>::parse_parameters (ParameterHandler &prm)
+    {
+      AssertThrow (dim==3, ExcMessage ("The 'radial linear with tidal potential' gravity model "
+                                       "can only be used in 3D."));
+      this->RadialLinear<dim>::parse_parameters (prm);
+
+      prm.enter_subsection("Gravity model");
+      {
+        prm.enter_subsection("Radial linear with tidal potnetial");
+        {
+          M_p = prm.get_double ("Mass of planet");
+          a_s = prm.get_double ("Semimajor axis of satellite");
+          b_NSR = prm.get_double ("Rate of nonsynchronous rotation");
+        }
+        prm.leave_subsection ();
+      }
+      prm.leave_subsection ();
+      // This effect of tidal potential only works if the geometry is derived from
+      // a spherical model (i.e. a sphere, spherical shell or chunk)
+      if (Plugins::plugin_type_matches<const GeometryModel::Sphere<dim>>(this->get_geometry_model()))
+        {
+          R1 = Plugins::get_plugin_as_type<const GeometryModel::Sphere<dim>>
+               (this->get_geometry_model()).radius();
+        }
+      else if (Plugins::plugin_type_matches<const GeometryModel::SphericalShell<dim>>(this->get_geometry_model()))
+        {
+          R1 = Plugins::get_plugin_as_type<const GeometryModel::SphericalShell<dim>>
+               (this->get_geometry_model()).outer_radius();
+        }
+      else if (Plugins::plugin_type_matches<const GeometryModel::Chunk<dim>>(this->get_geometry_model()))
+        {
+          R1 = Plugins::get_plugin_as_type<const GeometryModel::Chunk<dim>>
+               (this->get_geometry_model()).outer_radius();
+        }
+      else if (Plugins::plugin_type_matches<const GeometryModel::EllipsoidalChunk<dim>>(this->get_geometry_model()))
+        {
+          const auto &gm = Plugins::get_plugin_as_type<const GeometryModel::EllipsoidalChunk<dim>>
+                           (this->get_geometry_model());
+          // TODO
+          // If the eccentricity of the EllipsoidalChunk is non-zero, the radius can vary along a boundary,
+          // but the maximal depth is the same everywhere and we could calculate a representative pressure
+          // profile. However, it requires some extra logic with ellipsoidal
+          // coordinates, so for now we only allow eccentricity zero.
+          // Using the EllipsoidalChunk with eccentricity zero can still be useful,
+          // because the domain can be non-coordinate parallel.
+
+          AssertThrow(gm.get_eccentricity() == 0.0,
+                      ExcNotImplemented("This plugin cannot be used with a non-zero eccentricity. "));
+
+          R1 = gm.get_semi_major_axis_a();
+        }
+      else
+        {
+          Assert (false, ExcMessage ("This initial condition can only be used if the geometry "
+                                     "is a sphere, a spherical shell, a chunk or an "
+                                     "ellipsoidal chunk."));
+          R1 = numbers::signaling_nan<double>();
+        }
+    }
   }
 }
 
@@ -207,5 +365,17 @@ namespace aspect
                                   "is positive) and the magnitude changes linearly with depth. The "
                                   "magnitude of gravity at the surface and bottom is read from the "
                                   "input file in a section ``Gravity model/Radial linear''.")
+
+    ASPECT_REGISTER_GRAVITY_MODEL(RadialLinearWithTidalPotential,
+                                  "radial linear with tidal potential",
+                                  "A gravity model that is the sum of the `radial linear' model "
+                                  "(which is radial, pointing inward if the gravity "
+                                  "is positive, and a magnitude that changes linearly with depth), "
+                                  "and a term that results from a tidal gravity potential and that "
+                                  "leads to a gravity field that varies with latitude and longitude."
+                                  "The magnitude of gravity for the linear radial part is read from the "
+                                  "input file in a section `Gravity model/Radial linear'; the "
+                                  "parameters that describe the tidal potential contribution are read "
+                                  "from a section ``Gravity model/Radial linear with tidal potential''.")
   }
 }

tests/gravity_tidal_potential.prm

@@ -0,0 +1,100 @@
+# This parameter file tests the gravity model plugin for a case where the
+# tidal potential by flattening and non-synchronnous rotation changes gravity with position and time.
+# 
+# The equation implemented in this heating model is from Tobie et al. (2025) (https://doi.org/10.1007/s11214-025-01136-y),
+# which is defined as:
+# g = - magnitude - gradient ( - tidal potential ).
+# potential = 3 G M_p R_s^2 / 2 a_s^3 r^2 (Tstar + T0)
+# Tstar = 1/6 *(1-3*cos(theta)^2) and T0=1/2sin(theta)^2*cos(2*lambda + 2*b*t),
+# where G: gravitational constant, M_p: mass of planet, R_s: radius of satellite, a_s: semimajor axis of satellite's orbit, b = angular rate of nonsynchronous rotation.
+# r, theta and lambda are radial distance, polar angle and azimuthal angle, respectively.
+#
+# Model shows the Europa's icy shell without conduction in simpler model.
+# Visualization 'gravity' shows gravity distribution.
+
+set Dimension                              = 3
+set Use years instead of seconds           = true
+set End time                               = 1e4
+
+set Output directory                       = gravity_tidal_potential
+
+set Maximum first time step                    = 1e3
+set CFL number                                 = 0.8
+set Maximum time step                          = 1e3
+
+
+set Pressure normalization                 = surface
+set Surface pressure                       = 0
+
+
+subsection Geometry model
+  set Model name = spherical shell
+  subsection Spherical shell
+    set Outer radius = 1560800
+    set Inner radius = 1460800
+    set Opening angle = 360
+  end
+end
+
+
+subsection Initial temperature model
+  set Model name = function
+  subsection Function
+    set Coordinate system = spherical
+    set Variable names = r, phi,theta
+    set Function expression = 100
+  end
+end
+
+
+subsection Boundary velocity model
+  set Zero velocity boundary indicators = top, bottom
+end
+
+
+subsection Gravity model
+  set Model name = radial linear with tidal potential
+  subsection Radial linear
+    set Magnitude at surface = 1.3
+    set Magnitude at bottom = 1.3
+  end
+  subsection Radial linear with tidal potnetial
+  end
+end
+
+
+subsection Material model
+  set Model name = simpler
+  subsection Simpler model
+    set Reference density = 917
+    set Reference specific heat = 2110
+    set Reference temperature = 100
+    set Thermal conductivity = 0 #1.93
+    set Thermal expansion coefficient = 0 #1.6e-4
+    set Viscosity = 1e20
+  end
+end
+
+
+subsection Formulation
+  set Formulation = Boussinesq approximation
+end
+
+
+subsection Mesh refinement
+  set Initial global refinement                = 0
+  set Initial adaptive refinement              = 0
+  set Time steps between mesh refinement       = 0
+end
+
+
+subsection Postprocess
+  set List of postprocessors = velocity statistics, temperature statistics, visualization, basic statistics, \
+			       pressure statistics, material statistics
+
+  subsection Visualization
+    set Time between graphical output = 1e3
+    set Output format                 = vtu
+    set List of output variables      = material properties, strain rate, shear stress, stress, nonadiabatic pressure, gravity
+  end
+end