Skip to content

Commit 2cfde4d

Browse files
authored
Merge pull request geodynamics#6482 from hyunseong96/th_latitude2
Adding latitudinal variation to Tidal heating model
2 parents d25f9d8 + 348317e commit 2cfde4d

File tree

6 files changed

+402
-6
lines changed

6 files changed

+402
-6
lines changed
Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,5 @@
11
Added: ASPECT now has a new heating model plugin called 'tidal heating' for diurnal tides.
22
This plugin is useful for modeling long-term interior evolution in moons orbiting a large planet.
3+
Distribution of tidal strain rate can be selected between 'constant' and 'latitudinal variation'.
34
<br>
45
(Hyunseong Kim, 2025/06/13)

include/aspect/heating_model/tidal_heating.h

Lines changed: 28 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -23,6 +23,8 @@
2323
#define _aspect_heating_model_tidal_heating_h
2424

2525
#include <aspect/heating_model/interface.h>
26+
#include <aspect/simulator_access.h>
27+
2628

2729
namespace aspect
2830
{
@@ -34,8 +36,15 @@ namespace aspect
3436
* @ingroup HeatingModels
3537
*/
3638
template <int dim>
37-
class TidalHeating : public Interface<dim>
39+
class TidalHeating : public Interface<dim>, public ::aspect::SimulatorAccess<dim>
3840
{
41+
public:
42+
/**
43+
* Initialize function, which sets the start time and
44+
* start timestep of tidal heating.
45+
*/
46+
void initialize() override;
47+
3948
public:
4049
/**
4150
* Return the tidal heating terms.
@@ -74,10 +83,28 @@ namespace aspect
7483
* time-averaged strain rate = constant_tidal_strain_rate
7584
* tidal frequency = tidal_frequency
7685
* elastic shear modulus = elastic_shear_modulus
86+
*
87+
* strain_rate_distribution lets selection for distribution of tidal strain rate.
88+
* If 'constant', the tidal strain rate is fixed to 'Constant tidal strain rate'.
89+
* If 'latitudinal variation', 'Maximum tidal strain rate' and 'Minimum tidal strain rate' are used.
90+
* Then, the equation follows as (maximum_tidal_strain_rate - minimum_tidal_strain_rate)*cos(theta/2)+(maximum_tidal_strain_rate+minimum_tidal_strain_rate)/2.
91+
* This is shown in Fig.3 of Nimmo et al. (2007) (https://doi.org/10.1016/j.icarus.2007.04.021).
7792
*/
7893
double tidal_frequency;
7994
double elastic_shear_modulus;
8095
double constant_tidal_strain_rate;
96+
97+
/**
98+
* Specify the distribution of time-averaged tidal strain rate.
99+
*/
100+
enum StrainRateDistribution
101+
{
102+
constant,
103+
latitudinal_variation
104+
} strain_rate_distribution;
105+
106+
double maximum_tidal_strain_rate;
107+
double minimum_tidal_strain_rate;
81108
};
82109
}
83110
}

source/heating_model/tidal_heating.cc

Lines changed: 76 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -21,30 +21,73 @@
2121

2222
#include <aspect/heating_model/tidal_heating.h>
2323

24+
#include <aspect/simulator.h>
25+
#include <aspect/geometry_model/interface.h>
26+
#include <aspect/geometry_model/chunk.h>
27+
#include <aspect/geometry_model/two_merged_chunks.h>
28+
#include <aspect/geometry_model/ellipsoidal_chunk.h>
29+
#include <aspect/geometry_model/sphere.h>
30+
#include <aspect/geometry_model/spherical_shell.h>
2431

2532
namespace aspect
2633
{
2734
namespace HeatingModel
2835
{
36+
template <int dim>
37+
void
38+
TidalHeating<dim>::initialize ()
39+
{
40+
if (strain_rate_distribution == latitudinal_variation)
41+
{
42+
// Only spherical geometries will be accepted to have heating rate varying by latitude.
43+
AssertThrow((Plugins::plugin_type_matches<GeometryModel::Chunk<dim>>(this->get_geometry_model()) ||
44+
Plugins::plugin_type_matches<GeometryModel::TwoMergedChunks<dim>>(this->get_geometry_model()) ||
45+
Plugins::plugin_type_matches<GeometryModel::EllipsoidalChunk<dim>>(this->get_geometry_model()) ||
46+
Plugins::plugin_type_matches<GeometryModel::Sphere<dim>>(this->get_geometry_model()) ||
47+
Plugins::plugin_type_matches<GeometryModel::SphericalShell<dim>>(this->get_geometry_model())) &&
48+
dim==3,
49+
ExcMessage("Latitudinal variation of strain rate can only be used with 3-dimensional geometry models that "
50+
"have either a spherical or ellipsoidal natural coordinate system."));
51+
}
52+
}
53+
54+
55+
2956
template <int dim>
3057
void
3158
TidalHeating<dim>::
32-
evaluate (const MaterialModel::MaterialModelInputs<dim> &/*material_model_inputs*/,
59+
evaluate (const MaterialModel::MaterialModelInputs<dim> &material_model_inputs,
3360
const MaterialModel::MaterialModelOutputs<dim> &material_model_outputs,
3461
HeatingModel::HeatingModelOutputs &heating_model_outputs) const
3562
{
3663
/** *
3764
* H equation is from Tobie et al. (2003) (https://doi.org/10.1029/2003JE002099)
3865
* H= 2*(viscosity)*(time-averaged tidal strain rate)^2/(1+((viscosity)*(tidal frequency)/(elastic shear modulus))^2))
3966
* viscosity = material_model_outputs.viscosities
40-
* time-averaged strain rate = constant_tidal_strain_rate
67+
* time-averaged strain rate at certain location = local_strain_rate
4168
* tidal frequency = tidal_frequency
4269
* elastic shear modulus = elastic_shear_modulus
70+
* If 'Use latitudinal variation of strain rate' is true, local_strain_rate is calculated
71+
* with one period cosine function having maximum tidal strain rate at poles and minimum tidal strain rate at equator.
72+
* This variation of tidal strain rate follows as (maximum_tidal_strain_rate - minimum_tidal_strain_rate)*cos(theta/2)+(maximum_tidal_strain_rate+minimum_tidal_strain_rate)/2
73+
* where theta is polar angle from spherical coordinates.
74+
* If false, constant_tidal_strain_rate is used.
4375
*/
4476
for (unsigned int q=0; q<heating_model_outputs.heating_source_terms.size(); ++q)
4577
{
46-
heating_model_outputs.heating_source_terms[q] = 2 * material_model_outputs.viscosities[q] *constant_tidal_strain_rate * constant_tidal_strain_rate
47-
/ ( 1 + ( (tidal_frequency * material_model_outputs.viscosities[q]) / elastic_shear_modulus ) * ( (tidal_frequency * material_model_outputs.viscosities[q]) / elastic_shear_modulus ) );
78+
double local_tidal_strain_rate = 0.;
79+
if (strain_rate_distribution == latitudinal_variation)
80+
{
81+
const Point<dim> position = material_model_inputs.position[q];
82+
const double theta = Utilities::Coordinates::cartesian_to_spherical_coordinates(position)[dim-1];
83+
local_tidal_strain_rate = ( maximum_tidal_strain_rate - minimum_tidal_strain_rate ) / 2. * std::cos( theta * 2. ) + ( maximum_tidal_strain_rate + minimum_tidal_strain_rate ) / 2.;
84+
}
85+
else if (strain_rate_distribution == constant)
86+
{
87+
local_tidal_strain_rate = constant_tidal_strain_rate;
88+
}
89+
heating_model_outputs.heating_source_terms[q] = 2. * material_model_outputs.viscosities[q] * local_tidal_strain_rate * local_tidal_strain_rate
90+
/ ( 1. + ( (tidal_frequency * material_model_outputs.viscosities[q]) / elastic_shear_modulus ) * ( (tidal_frequency * material_model_outputs.viscosities[q]) / elastic_shear_modulus ) );
4891
heating_model_outputs.lhs_latent_heat_terms[q] = 0.0;
4992
}
5093
}
@@ -85,6 +128,23 @@ namespace aspect
85128
"Time-averaged strain rate by a lunar diurnal tide that is simplified to be constant regardless of location. "
86129
"Default value is for the diurnal tide of Europa from Tobie et al. (2003). "
87130
"Units: 1/s.");
131+
prm.declare_entry ("Custom distribution of tidal strain rate", "constant",
132+
Patterns::Selection("constant|latitudinal variation"),
133+
"Choose how the time-averaged tidal strain rate is distributed. "
134+
"If 'constant', the tidal strain rate is fixed to 'Constant tidal strain rate'. "
135+
"If 'latitudinal variation', 'Maximum tidal strain rate' and 'Minimum tidal strain rate' are used.");
136+
prm.declare_entry ("Maximum tidal strain rate", "2.81e-10",
137+
Patterns::Double (0),
138+
"Maximum time-averaged tidal strain rate by lunar diurnal tide at the poles. "
139+
"This parameter will be used when 'Custom distribution of tidal strain rate' is 'latitudinal variation'. "
140+
"Default value is for Europa at pole from Nimmo et al. (2007). "
141+
"Units: 1/s.");
142+
prm.declare_entry ("Minimum tidal strain rate", "1.67e-10",
143+
Patterns::Double (0),
144+
"Minimum time-averaged tidal strain rate by lunar diurnal tide at the equator. "
145+
"This parameter will be used when 'Custom distribution of tidal strain rate' is 'latitudinal variation'. "
146+
"Default value is for Europa at Equator from Nimmo et al. (2007). "
147+
"Units: 1/s.");
88148
}
89149
prm.leave_subsection();
90150
}
@@ -104,6 +164,14 @@ namespace aspect
104164
tidal_frequency = prm.get_double ("Tidal frequency");
105165
elastic_shear_modulus = prm.get_double ("Elastic shear modulus");
106166
constant_tidal_strain_rate = prm.get_double ("Constant tidal strain rate");
167+
if (prm.get ("Custom distribution of tidal strain rate") == "constant")
168+
strain_rate_distribution = constant;
169+
else if (prm.get ("Custom distribution of tidal strain rate") == "latitudinal variation")
170+
strain_rate_distribution = latitudinal_variation;
171+
else
172+
AssertThrow (false, ExcMessage ("Not a valid custom distribution of tidal strain rate."));
173+
maximum_tidal_strain_rate = prm.get_double ("Maximum tidal strain rate");
174+
minimum_tidal_strain_rate = prm.get_double ("Minimum tidal strain rate");
107175
}
108176
prm.leave_subsection();
109177
}
@@ -121,8 +189,11 @@ namespace aspect
121189
ASPECT_REGISTER_HEATING_MODEL(TidalHeating,
122190
"tidal heating",
123191
"A tidal heating implementation related to diurnal tides. "
124-
"The equation currently ignores regional (radial/lateral) changes. "
192+
"The default equation ignores regional (radial/lateral) changes. "
125193
"This equation is the Eq.12 from Tobie et al. (2003) (https://doi.org/10.1029/2003JE002099). "
194+
"Selecting 'latitudinal variation' from 'Custom distribution of tidal strain rate' allows simplified latitudinal variation "
195+
"with cosine function, 'Maximum tidal strain rate' and 'Minimum tidal strain rate. "
196+
"Latitudinal variation of tidal strain rate is shown in Fig.3 from Nimmo et al. (2007) (https://doi.org/10.1016/j.icarus.2007.04.021). "
126197
"Unit: W/m^3.")
127198
}
128199
}

tests/tidal_heating_latitude.prm

Lines changed: 112 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,112 @@
1+
# This parameter file tests the tidal heating plugin for a case where the
2+
# viscosity is constant and the tidal heating is dependent on latitude.
3+
#
4+
# The equation implemented in this heating model is from Tobie et al. (2003) (https://doi.org/10.1029/2003JE002099),
5+
# which is defined as:
6+
# H= 2*(viscosity)*(time-averaged tidal strain rate)^2/(1+((viscosity)*(tidal frequency)/(shear modulus))^2)), where
7+
# viscosity is the viscosity derived from the material model at every point (constant in this test - 1e14 Pa s)
8+
# The latitudinal variation of (time-averaged tidal strain rate) is simplified with cosine function between maximum tidal strain rate and minimum tidal strain rate.
9+
# The variation can be found in Fig.3 of Nimmo et al. (2007) (https://doi.org/10.1016/j.icarus.2007.04.021).
10+
#
11+
# The governing equation in the model is simplified as (density)*(specific heat capacity)*dT/dt=H, as unit of H is W/m^3.
12+
# As expected, temperature increases with time as the convective and conductive processes are not active.
13+
#
14+
# Analytical values are 2.85946709e+02 K and 1.65676318e+02 K at the pole and equator, respectively.
15+
# These values align with ASPECT's results within numerical accuracy. ASPECT's results are 2.86027243e+02 K and 1.65661824e+02 K at poles and equators, respectively.
16+
17+
set Dimension = 3
18+
set Use years in output instead of seconds = true
19+
set End time = 1e6
20+
21+
set Output directory = tidal_heating_latitude
22+
23+
set Maximum first time step = 1e5
24+
set CFL number = 0.8
25+
set Maximum time step = 1e5
26+
27+
28+
set Pressure normalization = surface
29+
set Surface pressure = 0
30+
31+
32+
subsection Geometry model
33+
set Model name = spherical shell
34+
subsection Spherical shell
35+
set Outer radius = 1560800
36+
set Inner radius = 1460800
37+
set Opening angle = 360
38+
end
39+
end
40+
41+
42+
subsection Initial temperature model
43+
set Model name = function
44+
subsection Function
45+
set Coordinate system = spherical
46+
set Variable names = r, phi,theta
47+
set Function expression = 100
48+
end
49+
end
50+
51+
52+
subsection Boundary velocity model
53+
set Zero velocity boundary indicators = top, bottom
54+
end
55+
56+
57+
subsection Gravity model
58+
set Model name = radial constant
59+
60+
subsection Radial constant
61+
set Magnitude = 0 #1.3
62+
end
63+
end
64+
65+
66+
subsection Material model
67+
set Model name = simpler
68+
subsection Simpler model
69+
set Reference density = 917
70+
set Reference specific heat = 2110
71+
set Reference temperature = 100
72+
set Thermal conductivity = 0 #1.93
73+
set Thermal expansion coefficient = 0 #1.6e-4
74+
set Viscosity = 1e14
75+
end
76+
end
77+
78+
79+
subsection Heating model
80+
set List of model names = tidal heating
81+
subsection Tidal heating
82+
set Tidal frequency = 2.048e-5
83+
set Elastic shear modulus = 3.3e9
84+
set Custom distribution of tidal strain rate = latitudinal variation
85+
set Maximum tidal strain rate = 2.81e-10
86+
set Minimum tidal strain rate = 1.67e-10
87+
end
88+
end
89+
90+
91+
subsection Formulation
92+
set Formulation = Boussinesq approximation
93+
end
94+
95+
96+
subsection Mesh refinement
97+
set Initial global refinement = 1
98+
set Initial adaptive refinement = 0
99+
set Time steps between mesh refinement = 0
100+
end
101+
102+
103+
subsection Postprocess
104+
set List of postprocessors = velocity statistics, temperature statistics, visualization, basic statistics, \
105+
pressure statistics, material statistics, heating statistics
106+
107+
subsection Visualization
108+
set Time between graphical output = 1e5
109+
set Output format = vtu
110+
set List of output variables = material properties, strain rate, shear stress, stress, nonadiabatic pressure, heating
111+
end
112+
end

0 commit comments

Comments
 (0)