Physics-Based Simulation

Semi-Implicit Discretization

However, challenges still remain on incorporating friction into the optimization time integration. A major problem is that friction is not a conservative force and there is no well-defined potential such that taking the opposite of its gradient produces the frictional force. In other words, implicit friction force is not integrable. Without a potential energy, backtracking line search could not be performed, and thus guarantees on the stability and convergence of the optimization will be broken.

In fact, whether a force has well-defined potential energy really depends on the temporal discretization. For example, with explicit time integration, any force $f$ is constant within a time step and it has a potential energy $-f^{T}x$. Taking this inspiration, we could make friction force integrable with a smarter temporal discretization. Making friction force constant within a time step would certainly restrict the size of the time step to obtain high quality results. Therefore, we discretize part of the friction force explicitly and formulate an integrable semi-implicit friction force.

Following IPC, we fix the normal force magnitude $\lambda$ (the ones only used in calculating friction) and the tangent operator $T$ during the nonlinear optimization to the value in the last time step $n$: ${\lambda }^n=\lambda (x^n)$, and $T^n=T(x^n)$, which then makes the friction force integrable with a potential energy $$P_f(x)=\sum _k\mu {\lambda }_k^n{f}_0(\Vert {\bar{\mathbf{v}}}_k\hat{h}\Vert ),\text{\hspace{1em}(9.2.1)}$$ where ${\bar{\mathbf{v}}}_k=(T_k^n{)}^{T}v$, $\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}\text{\hspace{1em}(9.2.2)}$$ so that $f_0^{\prime }(y)=f_1(y/\hat{h})$. Here $\hat{h}$ is a constant multiple of the time step size $h$ for most linear (multi-)step time integration methods including implicit Euler and higher-order backward difference formulas, etc. Then, taking the gradient of Equation (9.2.1) w.r.t. $x$ we obtain $$-\nabla P_f(x)=-\sum _k\mu {\lambda }_k^n{T}_k^n{f}_1(\Vert {\bar{\mathbf{v}}}_k\Vert )\mathbf{s}({\bar{\mathbf{v}}}_k),\text{\hspace{1em}(9.2.3)}$$ which is a semi-implicit discretization of our mollified friction force with explicit terms ${\lambda }_k^n$ and $T_k^n$. The Hessian of $P_f$ can be calculated as $$\begin{array}{rl}& {\nabla }^2{P}_f(x)\\ =& \sum _k\mu {\lambda }_k^n{T}_k^n(\frac{f_1^{\prime }(\Vert {\bar{\mathbf{v}}}_k\Vert )\Vert {\bar{\mathbf{v}}}_k\Vert -f_1(\Vert {\bar{\mathbf{v}}}_k\Vert )}{\Vert {\bar{\mathbf{v}}}_k{\Vert }^3}{\bar{\mathbf{v}}}_k{\bar{\mathbf{v}}}_k^T+\frac{f_1(\Vert {\bar{\mathbf{v}}}_k\Vert )}{\Vert {\bar{\mathbf{v}}}_k\Vert }{\mathbf{I}}_3){T_k^n}^T\frac{\partial v}{\partial x}.\end{array}\text{\hspace{1em}(9.2.4)}$$

Remark 9.2.1. In the friction gradient and Hessian expression (Equation (9.2.3) and Equation (9.2.4)), there are $\Vert {\mathbf{v}}_k\Vert$ in the denominators, which could be $0$ when there is no relative sliding motion at a contact point. To avoid division by $0$ during the computation, for friction gradient, we can derive $$\frac{f_1(\Vert {\bar{\mathbf{v}}}_k\Vert )}{\Vert {\bar{\mathbf{v}}}_k\Vert }=\{\begin{array}{ll}-\frac{\Vert {\bar{\mathbf{v}}}_k\Vert }{{ϵ}_v^2}+\frac{2}{{ϵ}_v},& \Vert {\bar{\mathbf{v}}}_k\Vert \in [0,{ϵ}_v)\\ 1/\Vert {\bar{\mathbf{v}}}_k\Vert ,& \Vert {\bar{\mathbf{v}}}_k\Vert \ge {ϵ}_v,\end{array}\text{\hspace{1em}(9.2.5)}$$ which is well-defined everywhere, and so we obtain $$-\nabla P_{f,k}(x)=-\mu {\lambda }_k^n{T}_k^n\frac{f_1(\Vert {\bar{\mathbf{v}}}_k\Vert )}{\Vert {\bar{\mathbf{v}}}_k\Vert }{\bar{\mathbf{v}}}_k=0\text{when}\Vert {\bar{\mathbf{v}}}_k\Vert =0.$$ For friction Hessian, we can derive $$\frac{f_1^{\prime }(\Vert {\bar{\mathbf{v}}}_k\Vert )\Vert {\bar{\mathbf{v}}}_k\Vert -f_1(\Vert {\bar{\mathbf{v}}}_k\Vert )}{\Vert {\bar{\mathbf{v}}}_k{\Vert }^2}=\{\begin{array}{ll}-1/{ϵ}_v^2,& \Vert {\bar{\mathbf{v}}}_k\Vert \in [0,{ϵ}_v)\\ -1/\Vert {\bar{\mathbf{v}}}_k{\Vert }^2,& \Vert {\bar{\mathbf{v}}}_k\Vert \ge {ϵ}_v,\end{array}\text{\hspace{1em}(9.2.6)}$$ which is also well-defined everywhere, and since ${\bar{\mathbf{v}}}_k{\bar{\mathbf{v}}}_k^T/\Vert {\bar{\mathbf{v}}}_k\Vert =0$ when $\Vert {\bar{\mathbf{v}}}_k\Vert =0$, we know that $${\nabla }^2{P}_{f,k}(x)=\mu {\lambda }_k^n{T}_k^n\left(\frac{f_1(\Vert {\bar{\mathbf{v}}}_k\Vert )}{\Vert {\bar{\mathbf{v}}}_k\Vert }{\mathbf{I}}_3\right){T_k^n}^T\frac{\partial v}{\partial x}\text{when}\Vert {\bar{\mathbf{v}}}_k\Vert =0.$$

Remark 9.2.2. The friction formulation in this lecture is introduced slightly differently from the original IPC [Li et al. 2020] in 2 places:

1. We directly use the relative sliding velocity ${\mathbf{v}}_k$ rather than the relative sliding displacement ${\mathbf{u}}_k=\hat{h}{\mathbf{v}}_k$ in IPC as the input to the mollifier $f_1()$, and so our $f_1()$ differs from that in the IPC on $\hat{h}$ in the denominators. When time integration rules other than implicit Euler is applied (so $x^{n+1}-x^n\ne \hat{h}{v}^{n+1}$), calling ${\mathbf{u}}_k$ the relative sliding displacement is inappropriate and may cause confusions.
2. We did not introduce a tangent basis to express relative sliding velocity in the tangent space, because this is not necessary in computing the friction energy, gradient, and Hessian.