pydrex.viscosity

  1"""> PyDRex: Dynamic viscosity based on various proposed constitutive equations.
  2
  3The contribution from dislocation creep can be calculated either using an exponential
  4law and a power-law equation for low and high homologous temperature respectively
  5(see [Demouchy et al. 2023](https://doi.org/10.2138/gselements.19.3.151)),
  6or using a unified temperature-insensitive creep law
  7(see [Garel et al. 2020](https://doi.org/10.1016/j.epsl.2020.116243)).
  8In the first case, it is necessary to choose a threshold homologous temperature¹ at which
  9the switch is made using the high-temperature power-law.
 10The referenced paper recommends a value of 0.7.
 11
 12¹Homologous temperature is defined as T/Tₘ where T is temperature and Tₘ is the melting
 13(solidus) temperature. The melting temperature of olivine is inversely proportional to
 14iron content and water content, and increases with pressure.
 15
 16"""
 17
 18import numpy as np
 19from scipy import constants
 20from scipy import special as sp
 21
 22import pydrex.minerals as _minerals
 23from pydrex.core import DefaultParams
 24
 25
 26def frictional_yielding(
 27    strain_rate_invariant_II,
 28    pressure,
 29    surface_cohesion,
 30    friction_coefficient,
 31    reference_density,
 32    # depth,
 33    # temperature,
 34):
 35    """Calculate viscosity contribution from semi-brittle yielding."""
 36    # TODO: Check: enumerator should come out as ~ 2e8.
 37    # TODO: Use depth and temperature thresholds? Check Bickert 2021 and Taufner 2021.
 38    τ_y = np.min(friction_coefficient + surface_cohesion * pressure, 1e9)
 39    return τ_y / (2 * strain_rate_invariant_II)
 40
 41
 42def matrix_dislocation_unified(
 43    strain_rate_invariant_II,
 44    temperature,
 45    pressure,
 46    dimension,
 47    coefficients=DefaultParams().disl_coefficients,
 48    method="tanh",
 49):
 50    """Calculate viscosity contribution from dislocation creep.
 51
 52    Calculate viscosity contribution from dislocation creep using a unified flow law
 53    parametrised by 3 polynomials (7 coefficients) as per
 54    [Garel et al. 2020](https://doi.org/10.1016/j.epsl.2020.116243).
 55
 56    """
 57    a0, b0, a1, b1, a2, b2, c2 = coefficients
 58    A0 = a0 + b0 * temperature
 59    A1 = a1 + b1 * temperature
 60    A2 = a2 + b2 * temperature + c2 * temperature**2
 61
 62    f = None
 63    match method:
 64        case "tanh":
 65            f = np.tanh
 66        case "erf":
 67            f = sp.erf
 68        case "algebraic":
 69
 70            def f(x):
 71                return x / np.sqrt(1 + x**2)
 72
 73    σ_diff = A0 * (1 + f(A1 * (np.log10(strain_rate_invariant_II) - A2)))
 74    return _adjust_for_geometry(strain_rate_invariant_II, dimension, σ_diff)
 75
 76
 77def diffusion(
 78    strain_rate_invariant_II,
 79    temperature,
 80    pressure,
 81    dimension,
 82    prefactors=DefaultParams().diff_prefactors,
 83    activation_energies=DefaultParams().diff_activation_energies,
 84    activation_volumes=DefaultParams().diff_activation_volumes,
 85    mode: str = "both",
 86):
 87    """Calculate viscosity contribution from diffusion creep.
 88
 89    Calculate viscosity contribution from diffusion in the
 90    `pydrex.core.DeformationRegime.matrix_diffusion` regime,
 91    using an exponential flow law.
 92
 93    - `strain_rate_invariant_II` — second invariant of the strain rate¹
 94    - `temperature` — temperature in Kelvin
 95    - `pressure` — pressure in GPa
 96    - `dimension` — number of coordinates describing a point in the domain
 97    - `prefactors` — prefactors for the exponential and power laws (s⁻¹)
 98    - `activation_energies` — see `pydrex.core.DefaultParams.disl_activation_energies`
 99    - `activation_volumes` — see `pydrex.core.DefaultParams.disl_activation_volumes`
100    - `mode` — which kind of diffusive regime is active (`"matrix"`, `"boundary"` or
101        `"both"`; `"sliding"` is an alias for `"both"`)
102
103    ¹The strain rate invariants can be calculated from the strain rate tensor using
104    `pydrex.tensors.invariants_second_order` if necessary.
105
106    """
107    # The pre-multipliers of α₁=14 and α₂=14π are from Kohlstedt & Hansen, 2015:
108    # <http://dx.doi.org/10.1016/B978-0-444-53802-4.00042-7> (see eq. 37).
109    # Technically, I think those values are only valid for the "both"|"sliding" mode,
110    # but I don't think there will be much chance to find values for other modes.
111    # In the rare case that only one of the two diffusive mechanisms is active and there
112    # is no sliding, I don't think that these 'geometric terms' will change a lot.
113
114    # TODO: Kohlstedt & Hansen 2015 also divide the prefactor by R*T again...
115    # This doesn't appear in Garel 2020 (eq. 10) nor in Chris' original models.
116    # For now we'll skip it but I have no idea which option is correct.
117    # Garel 2020 use R(T + δT) instead where δT is some sort of adiabatic correction,
118    # but once again I'll omit this here for now.
119    σ_diff_matrix = (
120        14
121        * prefactors[0] ** -1
122        * np.exp(
123            (activation_energies[0] + pressure * activation_volumes[0])
124            / (constants.gas_constant * temperature)
125        )
126    )
127    σ_diff_boundary = (
128        14
129        * np.pi
130        * prefactors[1] ** -1
131        * np.exp(
132            (activation_energies[1] + pressure * activation_volumes[1])
133            / (constants.gas_constant * temperature)
134        )
135    )
136    match mode:
137        case "matrix":
138            return _adjust_for_geometry(
139                strain_rate_invariant_II, dimension, σ_diff_matrix
140            )
141        case "boundary":
142            return _adjust_for_geometry(
143                strain_rate_invariant_II, dimension, σ_diff_boundary
144            )
145        case "both" | "sliding":
146            return _adjust_for_geometry(
147                strain_rate_invariant_II, dimension, σ_diff_matrix + σ_diff_boundary
148            )
149        case _:
150            raise ValueError(
151                f"diffusion creep mode must be one of 'matrix', 'boundary' or 'both', not '{mode}'"
152            )
153
154
155def matrix_dislocation_split(
156    strain_rate_invariant_II,
157    temperature,
158    pressure,
159    dimension,
160    prefactors=DefaultParams().disl_prefactors,
161    activation_energy=DefaultParams().disl_activation_energy,
162    activation_volume=DefaultParams().disl_activation_volume,
163    lowtemp_switch=DefaultParams().disl_lowtemp_switch,
164    σ_Peierls=DefaultParams().disl_Peierls_stress,
165    dry_olivine=True,
166):
167    """Calculate viscosity contribution from dislocation creep.
168
169    Calculate viscosity contribution from dislocation in the
170    `pydrex.core.DeformationRegime.matrix_dislocation` regime,
171    using an exponential law below the `lowtemp_switch` homologous temperature and a
172    power law above it.
173
174    See e.g. equations 2 and 3 in
175    [Demouchy et al. 2023](https://doi.org/10.2138/gselements.19.3.151).
176
177    - `strain_rate_invariant_II` — second invariant of the strain rate¹
178    - `temperature` — temperature in Kelvin
179    - `pressure` — pressure in GPa
180    - `dimension` — number of coordinates describing a point in the domain
181    - `prefactors` — prefactors for the exponential and power laws (s⁻¹)
182    - `activation_energy` — see `pydrex.core.DefaultParams.disl_activation_energy`
183    - `activation_volume` — see `pydrex.core.DefaultParams.disl_activation_volume`
184    - `lowtemp_switch` — see `pydrex.core.DefaultParams.disl_lowtemp_switch`
185    - `σ_Peierls` — see `pydrex.core.DefaultParams.disl_Peierls_stress`
186    - `dry_olivine` — use empyrical parametrisation for dry or hydrous olivine
187
188    ¹The strain rate invariants can be calculated from the strain rate tensor using
189    `pydrex.tensors.invariants_second_order` if necessary.
190
191    """
192    # Uniaxial strain rate correction.
193    strain_rate = 2 * np.sqrt(3) * strain_rate_invariant_II
194    homologous_temperature = temperature / _minerals.peridotite_solidus(pressure)
195
196    if homologous_temperature < lowtemp_switch:
197        σ_diff = σ_Peierls * np.power(
198            1
199            - np.power(
200                (np.log(strain_rate) * constants.gas_constant * temperature)
201                / (-activation_energy * prefactors[0]),
202                1.5,  # superscript q (supplement says 1 ≤ q ≤ 2)
203            ),
204            0.5,  # superscript p (supplement says 0 ≤ p ≤ 1)
205        )
206    else:
207        if dry_olivine:
208            n = 3.5
209        else:
210            n = 3
211
212        σ_diff = np.power(
213            strain_rate
214            / (
215                prefactors[1]
216                * np.exp(
217                    (-activation_energy + pressure * activation_volume)
218                    / (constants.gas_constant * temperature)
219                )
220            ),
221            (1 / n),
222        )
223
224    return _adjust_for_geometry(strain_rate_invariant_II, dimension, σ_diff)
225
226
227# This correction is discussed but not included in the Garel 2020 equations.
228# The idea is that we correct for the fact that σ_diff ≠ σ and
229# strain_rate_invariant_II ≠ strain_rate, where the deviation depends on the geometry.
230# The actual constitutive equations are usually calibrated against uniaxial compression
231# experiments or such, which is why they are actually defined in terms of σ_diff and
232# strain_rate_invariant_II rather than σ and strain_rate.
233def _adjust_for_geometry(strain_rate_invariant_II, dimension, σ_diff):
234    strain_rate = 2 * np.sqrt(3) * strain_rate_invariant_II
235    if dimension == 2:
236        return σ_diff / (2 * np.sqrt(3) * strain_rate)
237    elif dimension == 3:
238        return σ_diff / (3 * strain_rate)
239    else:
240        raise ValueError(f"dimension must be 2 or 3, not {dimension}")