Skip to content

Bathytropic Elasticity

The bathytropic elasticity model represents spatially heterogeneous isotropic elasticity whose stiffness is initialised from the integration-point depth. The depth dependence is prescribed through profiles of shear-wave velocity and density. Once the Young's modulus has been initialised, it remains constant at the material point and is independent of the subsequent stress, strain, void ratio, and current depth.

Contents

Overview

Bathytropic elasticity converts an initial spatial coordinate into a local elastic stiffness. At an integration point with material position \(\boldsymbol{X}\), the intended constitutive field is

\[ E=E(\boldsymbol{X}), \]

while Poisson's ratio \(\nu\) is constant. The stress increment is calculated using the local, isotropic stiffness tensor

\[ \Delta\boldsymbol{\sigma} = \mathbb{C}\big(E(\boldsymbol{X}),\nu\big):\Delta\boldsymbol{\varepsilon}. \]

The spatial variation of \(E\) does not evolve during the analysis. Consequently, the model differs fundamentally from barotropic elasticity: bathytropic elasticity does not update stiffness from the current pressure or void ratio.

Initial Depth Field

Let \(x_v\) denote the vertical coordinate of an integration point and let \(z^{ref}\) denote the reference elevation at which the depth is zero. The depth is

\[ d=\max\!\left(z^{ref}-x_v,0\right). \]

The implementation selects the vertical coordinate from the stress-vector size:

\[ x_v= \begin{cases} \texttt{coords(2)}, & \text{for } \texttt{ntens}=4,\\[2mm] \texttt{coords(3)}, & \text{for } \texttt{ntens}=6. \end{cases} \]

Thus, the vertical axis is the global \(y\)-axis in two-dimensional analyses and the global \(z\)-axis in three-dimensional analyses. Points above \(z^{ref}\) are assigned zero depth.

To make the depth functions dimensionally explicit, define the dimensionless depth

\[ \bar d=\frac{d}{d^{ref}}, \qquad d^{ref}=1\ \mathrm{m}. \]

The code expressions 1+z and log(1+z) correspond to \(1+\bar d\) and \(\ln(1+\bar d)\) when coordinates are given in metres.

Shear-Wave Velocity and Density Profiles

The shear-wave velocity profile is

\[ v_s(d)=v_s^0\left(1+\bar d\right)^m, \]

where \(v_s^0\) is the shear-wave velocity at zero depth and \(m\) controls its depth dependence.

The density profile used to calibrate the stiffness is

\[ \rho(d)=\rho^0+\epsilon_\rho\ln\left(1+\bar d\right), \]

where \(\rho^0\) is the density at zero depth and \(\epsilon_\rho\) controls its logarithmic change with depth.

For the commonly intended case \(m\ge0\) and \(\epsilon_\rho\ge0\), both \(v_s\) and \(\rho\) increase monotonically with depth. The constitutive law nevertheless only requires that the resulting density and elastic moduli remain positive over the depth range of the model.

Elastic Moduli

With \(v_s\) in m/s and \(\rho\) in kg/m³, the shear modulus in Pa is

\[ G_{\mathrm{Pa}}(d)=\rho(d)v_s(d)^2. \]

The implementation returns stress in kPa and therefore uses

\[ G(d)=\frac{\rho(d)v_s(d)^2}{1000}. \]

Young's modulus follows from the isotropic relation

\[ E(d)=2G(d)(1+\nu). \]

The Lamé parameter and bulk modulus are

\[ \lambda(d)=\frac{\nu E(d)}{(1+\nu)(1-2\nu)} \]

and

\[ K(d)=\frac{E(d)}{3(1-2\nu)}. \]

A positive-definite isotropic elastic stiffness requires

\[ E(d)>0 \qquad\text{and}\qquad -1<\nu<\frac{1}{2}. \]

Constitutive Relation and Strain Energy

After initialisation, \(E\) and \(\nu\) are constant at each integration point. In the small-strain setting, the local constitutive relation is

\[ \boldsymbol{\sigma} = \lambda(\boldsymbol{X})\,\operatorname{tr}(\boldsymbol{\varepsilon})\boldsymbol{I} +2G(\boldsymbol{X})\boldsymbol{\varepsilon}, \]

apart from any prescribed initial stress.

The corresponding strain-energy density is

\[ \Psi(\boldsymbol{\varepsilon};\boldsymbol{X}) = \frac{1}{2}\lambda(\boldsymbol{X}) \left[\operatorname{tr}(\boldsymbol{\varepsilon})\right]^2 +G(\boldsymbol{X})\boldsymbol{\varepsilon}:\boldsymbol{\varepsilon}. \]

It follows that

\[ \boldsymbol{\sigma}=\frac{\partial\Psi}{\partial\boldsymbol{\varepsilon}}. \]

Therefore, dependence of the elastic constants on the initial material position does not destroy hyperelasticity: under small-strain kinematics, the model is a spatially inhomogeneous hyperelastic material. The one-time depth initialisation introduces no evolving constitutive history after \(E\) has been stored.

Stress Update and Elastic Jacobian

The implementation updates the stress vector incrementally as

\[ \boldsymbol{\sigma}_{n+1} = \boldsymbol{\sigma}_n + \mathbf{D}^e(E,\nu)\Delta\boldsymbol{\varepsilon}. \]

Using engineering shear strains, the three-dimensional elastic matrix is

\[ \mathbf{D}^e= \begin{bmatrix} \lambda+2G & \lambda & \lambda & 0 & 0 & 0\\ \lambda & \lambda+2G & \lambda & 0 & 0 & 0\\ \lambda & \lambda & \lambda+2G & 0 & 0 & 0\\ 0 & 0 & 0 & G & 0 & 0\\ 0 & 0 & 0 & 0 & G & 0\\ 0 & 0 & 0 & 0 & 0 & G \end{bmatrix}. \]

For ntens = 4, only the first four rows and columns are used. Because the stored \(E\) is independent of the strain increment, the returned Jacobian

\[ \frac{\partial\Delta\boldsymbol{\sigma}} {\partial\Delta\boldsymbol{\varepsilon}} = \mathbf{D}^e \]

is the exact material Jacobian after initialisation.

Initialisation Algorithm

At every constitutive call, the implementation performs the following operations:

  1. If matvar(1) == 0, determine the vertical coordinate from coords, calculate \(d\), \(v_s(d)\), \(\rho(d)\), \(G(d)\), and \(E(d)\), and store \(E(d)\) in matvar(1).
  2. Otherwise, read \(E\) directly from matvar(1).
  3. Assemble the isotropic elastic matrix from \(E\) and \(\nu\).
  4. Update the stress with \(\Delta\boldsymbol{\sigma}=\mathbf{D}^e\Delta\boldsymbol{\varepsilon}\).

This algorithm freezes the modulus after the first nonzero value has been stored. It also has the following direct consequences:

  • any nonzero value supplied initially in matvar(1) is treated as the Young's modulus and bypasses the depth calculation;
  • a parameter combination that gives exactly \(E=0\) can never complete the sentinel-based initialisation;
  • the supplied implementation supports only stress vectors with ntens = 4 or ntens = 6.

Material Parameters

Index Symbol Unit Description
1 \(v_s^0\) m/s Reference shear-wave velocity at zero depth.
2 \(\rho^0\) kg/m³ Reference density at zero depth, used to derive \(G\).
3 \(\epsilon_\rho\) kg/m³ Logarithmic density-profile coefficient.
4 \(m\) - Shear-wave-velocity depth exponent.
5 \(\nu\) - Poisson's ratio.
6 \(z^{ref}\) m Reference elevation defining zero depth.

Internal State Variable

The implementation uses one state entry:

  1. matvar(1) — depth-initialised Young's modulus \(E\) in kPa.

The state variable is constant after initialisation. It must be retained across a restart if the original stiffness field is to remain unchanged. If no named output alias is registered, it is available as sdv1.

Scope and Limitations

Initial versus current coordinates

The intended material field is most naturally written as \(E(\boldsymbol{X})\), where \(\boldsymbol{X}\) is the reference position. The supplied routine, however, initialises from the coords array on the first constitutive call. In analyses without geometric nonlinearity, these are the original coordinates. With geometric nonlinearity, they are current coordinates. Initialisation must therefore occur before deformation if the reference-depth field is required. Material activation at a later stage produces a stiffness based on the coordinates at first activation.

Finite-strain interpretation

The position-dependent small-strain law admits the strain-energy density given above. Under geometric nonlinearity, the routine is nevertheless integrated through numgeo's general incremental Cauchy-stress framework. It should therefore not be interpreted as a fully finite-strain hyperelastic formulation derived from a deformation-gradient energy potential.

Density in wave-propagation analyses

The density \(\rho(d)\) in this constitutive model is an auxiliary quantity used to construct \(G(d)\). It does not update the density used in the mass matrix or body-force calculation. For a single-phase solid with mass density \(\rho_m\) expressed in t/m³, the numerical shear-wave speed is

\[ c_s(d)=\sqrt{\frac{G(d)}{\rho_m}}. \]

The specified profile \(c_s(d)=v_s(d)\) is recovered only when the mass density is assigned consistently, namely \(\rho_m(d)=\rho(d)/1000\) t/m³ under the unit convention used by the implementation.