Physics-Based Simulation

Discretization and Approximation

From Equation (18.4.2), the friction force per unit area is defined as

$${\mathbf{T}}_F(\mathbf{X},t)=-\mu \Vert {\mathbf{T}}_C(\mathbf{X},t)\Vert f(\Vert {\mathbf{V}}_F(\mathbf{X},t)\Vert )\mathbf{s}({\mathbf{V}}_F(\mathbf{X},t)),$$

where $\mu$ is the friction coefficient, ${\mathbf{T}}_C$ is the normal contact force, and ${\mathbf{V}}_F$ is the relative sliding velocity. Here $\mathbf{s}({\mathbf{V}}_F)=\frac{{\mathbf{V}}_F}{\Vert {\mathbf{V}}_F\Vert }$ when $\Vert {\mathbf{V}}_F\Vert >0$, while $\mathbf{s}({\mathbf{V}}_F)$ takes any unit vector orthogonal to the normal $\mathbf{N}(\mathbf{X},t)$ when $\Vert {\mathbf{V}}_F\Vert =0$. Additionally, the friction scaling function $f$ is also nonsmooth with respect to ${\mathbf{V}}_F$, as $f(\Vert {\mathbf{V}}_F\Vert )=1$ when $\Vert {\mathbf{V}}_F\Vert >0$, and $f(\Vert {\mathbf{V}}_F\Vert )\in [0,1]$ when $\Vert {\mathbf{V}}_F\Vert =0$.

It is important to note that without temporal discretization, there is no potential energy for friction. However, similar to Frictional Contact, once we discretize the normal force magnitude and the tangent operator to the last time step and smoothly approximate the friction scaling function $f$, the friction force at the $(n+1)$-th time step becomes integrable with respect to $\mathbf{x}$, and we obtain

$${\mathbf{T}}_F^{n+1}(\mathbf{X})\approx -\frac{\partial D^{n+1}(\mathbf{X})}{\partial {\mathbf{x}}^{n+1}(\mathbf{X})}=-\frac{\partial (\mu \Vert {\mathbf{T}}_C^n(\mathbf{X})\Vert f_0(\Vert {\bar{\mathbf{V}}}_F^{n+1}(\mathbf{X})\hat{h}\Vert ))}{\partial {\mathbf{x}}^{n+1}(\mathbf{X})}.$$

Here, ${\bar{\mathbf{V}}}_F^{n+1}(\mathbf{X})=(\mathbf{I}-{\mathbf{N}}^n(\mathbf{X}){\mathbf{N}}^n(\mathbf{X}{)}^T)({\mathbf{V}}^{n+1}(\mathbf{X})-{\mathbf{V}}^{n+1}({\mathbf{X}}_2))$ is the approximate relative sliding velocity, where ${\mathbf{N}}^n$ and ${\mathbf{X}}_2$ are the normal direction and the point in contact with $\mathbf{X}$ in the last time step, $\hat{h}I=(\partial v/\partial x{)}^{-1}$, and

$$f_0(y)=\{\begin{array}{ll}-\frac{y^3}{3{ϵ}_v^2{\hat{h}}^2}+\frac{y^2}{{ϵ}_v\hat{h}}+\frac{{ϵ}_v\hat{h}}{3},& y\in [0,{ϵ}_v\hat{h});\\ y,& y\ge {ϵ}_v\hat{h}.\end{array}$$

Therefore, considering self-contact, the approximate friction potential over the entire boundary can be written as

$${\int }_{{\Gamma }_C}\frac{1}{2}\mu \Vert {\mathbf{T}}_C^n(\mathbf{X})\Vert f_0(\Vert {\bar{\mathbf{V}}}_F^{n+1}(\mathbf{X})\hat{h}\Vert )ds(\mathbf{X}),$$

where the $\frac{1}{2}$ scaling comes from double counting the friction between each pair of contact points in the integral (similar to the normal contact forces in Boundary Conditions and Frictional Contact).

After discretizing the boundary curves as polylines and approximating the max operator in the normal contact force component using summations (Piecewise Linear Boundaries), we similarly obtain the spatially discretized friction potential:

$${\int }_{{\Gamma }_C}\sum _{e\in \mathcal{E}-I(\mathbf{X})}\frac{1}{2}\mu \left(-\frac{\partial b(d^{\text{PE}}({\mathbf{x}}^n(\mathbf{X}),e),\hat{d})}{\partial d}\right)f_0(\Vert {\bar{\mathbf{V}}}_F^{n+1}(\mathbf{X},e)\hat{h}\Vert )ds(\mathbf{X}).$$

Here, $d^{\text{PE}}({\mathbf{x}}^n(\mathbf{X}),e)$ is the point-edge distance between ${\mathbf{x}}^n(\mathbf{X})$ and edge $e$ in the last time step, and ${\bar{\mathbf{V}}}_F^{n+1}(\mathbf{X},e)$ is the approximate relative sliding velocity of the point-edge pair with contact normal and the closest point discretized to the last time step (see next section for details).

If we choose boundary nodes as quadrature points to approximate the integral, we finally obtain our discrete friction potential:

$$\begin{array}{rl}{P}_f(x)& =\sum _{\hat{a}}{A}_{\hat{a}}\sum _{e\in \mathcal{E}-I({\mathbf{X}}_{\hat{a}})}\frac{1}{2}\mu \left(-\frac{\partial b(d^{\text{PE}}({\mathbf{x}}_{\hat{a}}^n,e),\hat{d})}{\partial d}\right)f_0(\Vert {\bar{\mathbf{V}}}_F^{n+1}({\mathbf{X}}_{\hat{a}},e)\hat{h}\Vert )\\ & =\sum _{k\in \{(\hat{a},e)\}}\mu {\lambda }_k^n{f}_0(\Vert {\bar{\mathbf{v}}}_k\hat{h}\Vert )\end{array}$$

where $A_{\hat{a}}=\frac{\Vert {\mathbf{X}}_{\hat{a}}-{\mathbf{X}}_{\hat{a}-1}\Vert +\Vert {\mathbf{X}}_{\hat{a}}-{\mathbf{X}}_{\hat{a}+1}\Vert }{2}$ is the integration weight. Denoting ${\bar{\mathbf{v}}}_k={\bar{\mathbf{V}}}_F^{n+1}({\mathbf{X}}_{\hat{a}},e)$ and ${\lambda }_k^n=\frac{1}{2}{A}_{\hat{a}}\left(-\frac{\partial b(d^{\text{PE}}({\mathbf{x}}_{\hat{a}}^n,e),\hat{d})}{\partial d}\right)$, the expression of $P_f$ agrees with the discrete form of Equation (9.2.1) we directly derived, except that here $k$ traverses all non-incident point-edge pairs on the boundary.

Based on this discrete form of the smoothed semi-implicit friction potential, we now need to determine how to calculate $\lambda$ and $\bar{\mathbf{v}}$ for point-edge pairs, implement the computation of the value, gradient, and Hessian of $P_f(x)$, and then incorporate them into the optimization.