Local contact coordinates

Contact geometry is given by local coordinates (situated at the contact point and having the \(x\)-axis along the contact normal) and also is responsible for computing relative velocity of contacting particles in local coordinates. This functionality is implemented by functors producing woo.dem.L6Geom (woo.dem.Cg2_Sphere_Sphere_L6Geom and others).

Local coordinates are defined by position (local origin) \(c\) (contact point) and rotation (local orientation) given as rotation (=orthonormal) matrix \(\mat{T}\). The first row of \(\mat{T}\) i.e. the local \(x\)-axis, is noted \(\vec{n}\) for brevity.

There is 6 generalized displacements (degrees of freedom) at the contact point:

relative displacements \(\vec{u}\)

the first (\(x\)) component is the normal displacement \(u_n\), the other two are tangential displacement vector \(\vec{u_t}\);

relative rotations \(\vec{\phi}\)

the \(x\)-component is the twist rotation \(\phi_x\), \(\phi_y\) and \(\phi_z\) are bending (in-plane) rotation along respective local axes.

All displacements define their respective rates noted \(\dot{\vec{u}}\) and \(\dot{\vec{\phi}}\).

Cg2_* functors are responsible for updating these values on every contact. Most of the computation is done incrementally.

Note

It must be carefully observed for which time-point is any value valid; this is show by using time-indices \(\prev{x}\equiv x^{(t-\Dt)}\), \(\pprev{x}\equiv x^{(t-\Dt/2)}\), \(\curr{x}\equiv x^{(t)}\), \(\nnext{x} \equiv x^{(t+\Dt/2)}\), \(\next{x} \equiv x^{(t+\Dt)}\). Because of the leap-frog integration scheme, even derivatives of position (position, acceleration, forces, torques) are known at full-steps while odd derivatives (velocity) are known at mid-steps.

Generic contact routine

Contact of any two shapes is handled by the same routine which supposes two spheres with radii \(r_i\), positions \(\vec{x}_i\), velocities \(\vec{v}_i\), angular velocities \(\vec{\omega}_i\). If the shapes in contact are not really spheres, we use “equivalent” sphere geometry which captures the configuration in question.

At every step, every contact directly computes (incomplete) current local coordinates: the contact point \(\curr{\vec{c}}\) and contact normal \(\curr{\vec{n}}\). The rest is computed indirectly. This computation is different for different combinations of shapes (sphere with sphere, wall, facet and so on).

Contact formation

When a contact is initially formed, current values of \(\curr{\vec{c}}\) and \(\curr{\vec{n}}\) are computed directly. Local axes (columns of \(\mat{T}\)) are initially defined as follows:

• local \(x\)-axis is contact normal \(\vec{n}\);

• local \(y\)-axis is positioned arbitrarily, but in a deterministic manner: aligned with the \(xz\) plane (if \(\vec{n}_y<\vec{n}_z\)) or \(xy\) plane (otherwise);

• local \(z\)-axis is perpendicular to both other axis \(\vec{z}_l=\vec{x}_l\times\vec{y}_l\).

Relative velocities are zero \(\dot{\curr{\vec{u}}}=\dot{\curr{\vec{\phi}}}=\vec{0}\), displacements \(\vec{u}\) and \(\vec{\phi}\) are also intially zero (normal displacement is often compute in a non-incremental way and may have non-zero initial value).

Todo

Is it ok that relative velocities are initially zero, and not updated until in the next step? Those could be computed directly, although with less precision, when the contact is formed.

Note

Quasi-constant values of \(\vec{u}_0\) (and \(\vec{\phi}_0\)) are stored as shifted origins of \(\vec{u}\) (and \(\vec{\phi}\)); therefore, the current value of displacement is always \(\curr{\vec{u}}-\vec{u}_0\) This is useful for easy implementation of plastic slipping in the normal/twist sense where the value is often computed non-incrementally.

Contact update

For a contact which already exists, the value of \(\curr{\vec{c}}\) and \(\curr{\vec{n}}\) is computed directly again. The contact is then updated to keep track of rigid motion of the contact (one that does not change mutual configuration of spheres) and mutual configuration changes. Some computations can be more or less approximate (trading performance for accuracy), which is controlled by approximation flags.

Rigid motion transforms local coordinate system and can be decomposed in rigid translation (affecting \(\vec{c}\)), and rigid rotation (affecting \(\mat{T}\)), which can be split in bending rotation \(\vec{o}_r\) perpendicular to the normal and twisting rotation \(\vec{o}_t\) parallel with the normal:

\[\pprev{\vec{o}_r}=\prev{\vec{n}}\times\curr{\vec{n}}.\]

Since velocities are known at previous midstep (\(t-\Dt/2\)), we consider mid-step normal

\[\begin{split}\pprev{\vec{n}}\begin{cases}=\frac{\prev{\vec{n}}+\curr{\vec{n}}}{2} & \text{(accurate solution)} \\ \approx\prev{\vec{n}} & \text{(with approximation flag set)}\end{cases}.\end{split}\]

For the sake of numerical stability, \(\pprev{\vec{n}}\) is re-normalized after being computed (unless prohibited by approximation flags).

Rigid rotation parallel with the normal is

\[\pprev{\vec{o}_t}=\pprev{\vec{n}}\left(\pprev{\vec{n}}\cdot\frac{\pprev{\vec{\omega}}_1+\pprev{\vec{\omega}}_2}{2}\right)\Dt.\]

Branch vectors \(\vec{b}_1\), \(\vec{b}_2\) (connecting \(\curr{\vec{x}}_1\), \(\curr{\vec{x}}_2\) with \(\curr{\vec{c}}\) are computed depending on noRatch (see details in Yade docs):

\begin{align*} \vec{b}_1&=\begin{cases} r_1 \curr{\vec{n}} & \mbox{with noRatch} \\ \curr{\vec{c}}-\curr{\vec{x}}_1 & \mbox{otherwise} \end{cases} \\ \vec{b}_2&=\begin{cases} -r_2\curr{\vec{n}} & \mbox{with noRatch} \\ \curr{\vec{c}}-\curr{\vec{x}}_2 & \mbox{otherwise} \end{cases} \\ \end{align*}

Relative velocity at \(\curr{\vec{c}}\) can be computed as

\[\pprev{\vec{v}_r}=(\pprev{\vec{\tilde{v}}_2}+\vec{\omega}_2\times\vec{b}_2)-(\vec{v}_1+\vec{\omega}_1\times\vec{b}_1)\]

where \(\vec{\tilde{v}}_2\) is \(\vec{v}_2\) without mean-field velocity gradient in periodic boundary conditions (see woo.core.Cell.homoDeform). In the numerial implementation, the normal part of incident velocity is removed (since it is computed directly) and replaced with with \(\pprev{\vec{v}_{r2}}=\pprev{\vec{v}_r}-(\pprev{\vec{n}}\cdot\pprev{\vec{v}_r})\pprev{\vec{n}}\).

Any vector \(\vec{a}\) expressed in global coordinates transforms during one timestep as

\[\curr{\vec{a}}=\prev{\vec{a}}+\pprev{\vec{v}_r}\Dt-\prev{\vec{a}}\times\pprev{\vec{o}_r}-\prev{\vec{a}}\times{\pprev{\vec{t}_r}}\]

where the increments have the meaning of relative shear, rigid rotation normal to \(\vec{n}\) and rigid rotation parallel with \(\vec{n}\). Local coordinate system orientation, rotation matrix \(\mat{T}\), is updated by columns (written as transpose for space reasons), i.e.

\[\begin{split}\curr{\mat{T}}=\begin{pmatrix} \curr{\vec{n}_x}, \curr{\vec{n}_y}, \curr{\vec{n}_z} \\ \prev{\mat{T}_{1,\bullet}}-\prev{\mat{T}_{1,\bullet}}\times\pprev{\vec{o}_r}-\prev{\mat{T}_{1,\bullet}}\times\pprev{\vec{o}_t} \\ \prev{\mat{T}_{2,\bullet}}-\prev{\mat{T}_{2,\bullet}}\times\pprev{\vec{o}_r}-\prev{\mat{T}_{,\bullet}}\times\pprev{\vec{o}_t} \end{pmatrix}^T\end{split}\]

This matrix is re-normalized (unless prevented approximation flags) and mid-step transformation is computed using quaternion spherical interpolation as

\[\begin{split}\pprev{\mat{T}}\begin{cases}=\mathrm{Slerp}\,\left(\prev{\mat{T}};\curr{\mat{T}};t=1/2\right) & \text{(accurate solution)} \\ \approx\prev{\mat{T}} & \text{(with approximation flag set)}\end{cases}.\end{split}\]

Finally, current generalized displacements in local coordinates are evaluated as

\begin{align*} \curr{\vec{u}}&=\prev{u}+{\pprev{\mat{T}}}^T\pprev{\vec{v}_r}\Dt, \\ \curr{\vec{\phi}}&=\prev{\vec{\phi}}+{\pprev{\mat{T}}}^T\Dt(\vec{\omega}_2-\vec{\omega}_1) \end{align*}

For the normal component, non-incremental evaluation is preferred if possible; for two spheres, this reads

\[\curr{\vec{u}_x}=|\curr{\vec{x}_2}-\curr{\vec{x}_1}|-(r_1+r_2).\]

Tip

Report issues or inclarities to github.