digger.analyze.cm_tools

Module Contents

Functions

CM_vel_acc

Calculate the center of mass velocity using spatial gradients.

API

digger.analyze.cm_tools.CM_vel_acc(x: numpy.ndarray, y: numpy.ndarray, H: numpy.ndarray, HU: numpy.ndarray, HV: numpy.ndarray, HM: numpy.ndarray, P: numpy.ndarray, ETA: numpy.ndarray, rhos: float, rhof: float, g: float = 9.81, drytolerance: float = 0.001, phi: float = 38.0, component: str = 'solid', aoi: None | numpy.ndarray = None, solidfrac: None | numpy.ndarray = None)

Calculate the center of mass velocity using spatial gradients.

Theory of center of mass velocity an acceleration estimation

Determining the motion of the center of mass of a landslide body is useful for energetic analysis or other applications, for example, landslide seismology.

The center of mass \(\vec{X}_{cm}(t)\) of a variably dense mixture of thickness \(h(x,y,t)\), occupying the region between the fixed bed \(b(x,y)\) and upper free-surface \(\eta(x,y,t)\) is

\[\begin{align}\label{CM_cont} \vec{X}_{cm}(t) = \frac{1}{M}\int_{\Omega}\rho(\vec{x},t)\vec{x}\,dV = \frac{1}{M}\int_{\Omega_\bot} \left[ \rho h x, \rho h y, \rho h\bar{z}\right]\,dx\,dy, \end{align}\]

where \(\rho(\vec{x},t)=\rho(x,y,t)\) is the bulk density (assumed to be independent of \(z\)), \(M\) is the total mass, and \(\bar{z} := \tfrac{1}{2h}\int^{\eta}_{b} z\,dz = \tfrac{1}{2}(\eta+b)= b+\tfrac{1}{2}h\).

Considering the mass \(M\) and domain \(\Omega\) to be fixed, the center of mass velocity, \(\vec{X}^\prime_{cm}(t)\), is given by

\[\begin{split}\begin{multline}\label{CMV_cont} \vec{X}^{\prime}_{cm}(t) = \frac{1}{M}\int_{\Omega_\bot} \left[ x(\rho h)_t , y(\rho h)_t , (\bar{z}\rho h)_t\right]\,dx\,dy\\ = \frac{1}{M}\int_{\Omega_\bot} \left[ x(\rho h)_t, y(\rho h)_t, \bar{z}(\rho h)_t+ \bar{z}_t\rho h \right]\,dx\,dy,\\ = \frac{1}{M}\int_{\Omega_\bot} \left[ x(\rho h)_t, y(\rho h)_t, (b+\tfrac{1}{2}h)(\rho h)_t+ \tfrac{1}{2}h_t\rho h \right]\,dx\,dy. \end{multline}\end{split}\]

Similarly, the center of mass acceleration, \(\vec{X}_{cm}^{\prime\prime}(t)\), is given by

\[\begin{split}\begin{multline}\label{CMA_cont} \vec{X}^{\prime\prime}_{cm}(t) = \frac{1}{M}\int_{\Omega_\bot} \left[ (\rho h)_{tt} x, (\rho h)_{tt} y, (\rho h\bar{z})_{tt}\right]\,dx\,dy,\\ = \frac{1}{M}\int_{\Omega_\bot} \left[ x(\rho h)_{tt}, y(\rho h)_{tt},\left( (b+\tfrac{1}{2}h)(\rho h)_t+ \tfrac{1}{2}h_t\rho h \right)_t\right]\,dx\,dy,\\ = \frac{1}{M}\int_{\Omega_\bot} \left[ x(\rho h)_{tt}, y(\rho h)_{tt},(b+\tfrac{1}{2}h)(\rho h)_{tt}+ \tfrac{1}{2}h_{tt}\rho h + h_t(\rho h)_{t} \right]\,dx\,dy. \end{multline}\end{split}\]

For a piecewise-constant function (\(h_{ij}(t), \rho_{ij}(t))\) on a uniform grid, with \(x_{ij} = x_0 + (i-\tfrac{1}{2})\Delta x\), and \(y_{ij} = y_0 + (j-\tfrac{1}{2})\Delta y\),

\[\begin{align}\label{CM} \vec{X}_{cm}(t)= \frac{\Delta x\Delta y}{M}\sum_{i,j} \left[ (x\rho h)_{ij}, (y\rho h)_{ij}, (\bar{z}\rho h)_{ij} \right], \end{align}\]

where \(M=\sum_{i,j} (\rho h)_{ij}\Delta x \Delta y\). Similarly, \(\vec{X}^{\prime}_{cm}(t)\) and \(\vec{X}^{\prime\prime}_{cm}(t)\) can be written as grid summations of the time derivatives of \((\rho h)_{ij}$ and $\bar{z}_{ij}\), as in the integrals above.

For a numerical solution the variables \(\rho\) and \(h\) are only calculated at discrete times. Calculating the center of mass velocity, \(\vec{X}^{\prime}_{cm}(t)\), or acceleration \(\vec{X}^{\prime\prime}_{cm}(t)\), by differencing the - discrete trajectory \(\vec{X}_{cm}(t)\) at times \(t=t^n, n=1,\dots,N\), is numerically problematic because error or noise in the data \(\vec{X}_{cm}(t)\) is magnified with each derivative. An alternative approach, the one used here, is to determine \(\vec{X}^{\prime}_{cm}(t^n)\) and \(\vec{X}^{\prime\prime}_{cm}(t^n)\) by using the numerical solution only at \(t^n\) together with the PDEs that govern the numerical solution. In other words, the time derivatives of \((\rho h)_{ij}\) are replaced by discrete spatial operators derived from the PDEs and applied to \((\rho h)_{ij}^n\).

Consider the following PDEs for a variably dense but incompressible mixture,

\[\begin{split}\begin{align}\label{SWE} h_t + (hu)_x + (hv)_y & =0,\\ (\rho h)_t + (\rho h u)_x + (\rho h v)_y &= 0, \\ (\rho hu)_t + \left( \tfrac{1}{2} \rho g h^2 + \rho hu^2\right)_x + (\rho h u v)_y &= -\rho ghb_x -\mathrm{sgn}(u)\tau^x,\\ (\rho hv)_t + (\rho h u v)_x + \left( \tfrac{1}{2} \rho g h^2 + \rho hv^2\right)_y &= -\rho ghb_y -\mathrm{sgn}(v)\tau^y, \end{align}\end{split}\]

where \(\tau^x\) and \(\tau^y\) are the shear resistance in the \(x\) and \(y\) directions respectively (we assume the solution is twice differentiable). From the first two equations,

\[\begin{split}\begin{align} h_t = -(h u)_x -( h v)_y,\quad\mathrm{and},\\ (\rho h)_t = -(\rho h u)_x -(\rho h v)_y. \end{align}\end{split}\]

The second time derivatives are replaced by utilizing the momentum equations,

\[\begin{split}\begin{multline}\label{hrho_tt} (h\rho)_{tt} = -((\rho h u)_{x}+(\rho h v)_{y})_t= -((\rho h u)_{xt}+(\rho h v)_{yt}) = -\left(((\rho h u)_{t})_x+((\rho h v)_{t})_y\right)=\\ \left(\left( \tfrac{1}{2} \rho g h^2 + \rho hu^2\right)_x + (\rho h u v)_y +\rho ghb_x +\mathrm{sgn}(u)\tau^x\right)_x\\ +\left((\rho h u v)_x + \left( \tfrac{1}{2} \rho g h^2 + \rho hv^2\right)_y +\rho ghb_y +\mathrm{sgn}(v)\tau^y \right)_y=\\ \left( \tfrac{1}{2} \rho g h^2 + \rho hu^2\right)_{xx} + (\rho h u v)_{yx} +\rho ghb_{xx}+ g(\rho h)_x b_{x} + (\mathrm{sgn}(u)\tau^x)_x\\ +(\rho h u v)_{xy} + \left( \tfrac{1}{2} \rho g h^2 + \rho hv^2\right)_{yy} +\rho ghb_{yy} +g(\rho h)_yb_y+(\mathrm{sgn}(v)\tau^y)_y. \end{multline}\end{split}\]

For the term \(h_{tt}\), we neglect local gradients in \(\rho\) and use the standard shallow-flow momentum equations, leading to:

\[\begin{split}\begin{multline}\label{h_tt} h_{tt} = \left( \tfrac{1}{2} g h^2 + hu^2\right)_{xx} + ( h u v)_{yx} + ghb_{xx}+ gh_x b_{x} + (\mathrm{sgn}(u)\tau^x/\rho)_x\\ +( h u v)_{xy} + \left( \tfrac{1}{2} g h^2 + hv^2\right)_{yy} +ghb_{yy} +g h_yb_y+(\mathrm{sgn}(v)\tau^y/\rho)_y. \end{multline}\end{split}\]

To summarize, we numerically evaluate the center of mass velocity and acceleration using the numerical solution at \(t^n\):

\[\begin{align}\label{CMV_sum} \vec{X}^{\prime}_{cm}(t^n) = \frac{\Delta x\Delta y}{M}\sum_{i,j} \left[ x(\rho h)_t, y(\rho h)_t, (b+\tfrac{1}{2}h)(\rho h)_t+ \tfrac{1}{2}h_t\rho h\right]_{ij}, \end{align}\]

and

\[\begin{align}\label{CMA_sum} \vec{X}^{\prime\prime}_{cm}(t^n) = \frac{\Delta x\Delta y}{M}\sum_{i,j} \left[ x(\rho h)_{tt}, y(\rho h)_{tt},(b+\tfrac{1}{2}h)(\rho h)_{tt}+ \tfrac{1}{2}h_{tt}\rho h + h_t(\rho h)_{t} \right]_{ij}, \end{align}\]

where the time derivatives appearing in the sums are approximated by spatial discretizations using the above derivations.

All units assumed to be meters-kilograms-seconds.

Inputs:

xnp.ndarray

1-d numpy array of the x-coordinates of grid cell centers. Algorithm assumes that x is increasing.

ynp.ndarray

1-d numpy array of the y-coordinates of grid cell centers. Algorithm assumes that y is increasing.

Hnp.ndarray

2-d numpy array of flow depth.

HUnp.ndarray

2-d numpy array of depth times x-directed velocity.

HV: np.ndarray

2-d numpy array of depth times y-directed velocity.

HM: np.ndarray

2-d numpy array of depth times solid volume fraction.

P: np.ndarray

2-d numpy array of basal pressure.

ETA: np.ndarray

2-d numpy array of the basal surface plus the flow depth

rhos: float

Solid-fraction density.

rhof: float

Fluid fraction density.

componentstr

Whether to consider all material (‘all’), the fluid fraction (‘fluid’), or the solid fraction (‘solid’).

g: float = 9.81

Gravitational acceleration.

drytolerance: float = 0.001

Dry tolerance for analysis. Only cells with depth greater than drytolerance are considered.

aoinp.ndarray

2-d numpy array of boolean values indicating an area of interest (True) and the rest of the area (False). If not provided all cells are considered.

solidfracnp.ndarray

2-d numpy array of boolean values indicating an area where the solid fraction is of interest. If not provided all cells are considered.

Outputs:

cmtuple of three floats

Size three tuple containing the x, y, and z coordinates of the center of mass.

cmvtuple of three floats

Size three tuple containing the x-directed, y-directed, and z-directed velocities of the center of mass.

cmatuple of three floats

Size three tuple containing the x-directed, y-directed, and z-directed accelerations of the center of mass.