Quasistatic Deformable Simulation

Continuum Mechanics

What is the continuous problem we want to solve before we think about discretization schemes?

\begin{aligned} \text{find}\quad & v(x) \\ \text{s.t.}\quad & \int_{\Omega_+} \nabla\cdot \sigma(v) dV + \int_{\Omega_+} \rho g dV = \int_{\partial \Omega_+} \lambda \cdot \hat{n}(v) dA \end{aligned}

\Omega=\{x | x \in\Omega\subset\mathbb{R}^3\}

\Omega_+=\{x_+ | x_+= x + v(x),\; x\in\Omega\}

Current configuration

Next configuration

Find the displacement field subject to quasistatic force balance

Elastic deformation

Body forces (gravity)

Boundary traction forces

Continuum Mechanics

\begin{aligned} \text{find}\quad & v(x) \\ \text{s.t.}\quad & \int_{\Omega_+} \nabla\cdot \sigma(v) dV + \int_{\Omega_+} \rho g dV = \int_{\partial \Omega_+} \lambda \cdot \hat{n}(v) dA \end{aligned}

In linear elasticity, the strain dependent forces are linearly dependent on strain, which is the gradient of our decision variable.

This equation can be solved with Galerkin's method (e.g. FEM, MPM)
But the boundary traction forces are tough to handle. Folks have used penalty-based methods to explicitly define contact forces.
But our application is a bit special since it involves large changes in boundary conditions. Can we take this formulation to what we had in the rigid-body case?

\begin{aligned} \sigma_{ij} & = K_{ijkl} \varepsilon_{kl} \qquad \varepsilon_{kl} & = \frac{\partial v_k}{\partial x_l} \qquad (\nabla\cdot\sigma)_i & = \frac{\partial \sigma_{ij}}{\partial x_j} = \frac{\partial}{\partial x_j} K_{ijkl}\frac{\partial v_k}{\partial x_l} \end{aligned}

Contribution 1. Energy optimization

\begin{aligned} \text{find}\quad & v(x) \\ \text{s.t.}\quad & \int_{\Omega_+} \nabla\cdot \sigma(v) dV + \int_{\Omega_+} \rho g dV = \int_{\partial \Omega_+} \lambda \cdot \hat{n}(v) dA \end{aligned}

We turn this into an optimization problem whose optimality conditions imply force balance.

\begin{aligned} \text{min}_{v(x)}\quad & \int_{\Omega_+}\Psi(v)dV + \int_{\Omega_+}\rho g v dV \\ \text{s.t.}\quad & \phi_{\Omega_+}(x+v) \geq 0 \quad \forall x\in\Omega \end{aligned}

\begin{aligned} \Psi(v) & = \frac{1}{2}K_{ijkl} \varepsilon_{ij}\varepsilon_{kl} \\ & = \frac{1}{2}K_{ikjl} \frac{\partial v_i}{\partial x_j}\frac{\partial v_k}{\partial x_l} \end{aligned}

Potential strain energy

Multiplier term

\begin{aligned} \delta\Psi(v) & = \frac{1}{2}K_{ijkl}\frac{\delta v_i}{\delta x_j}\frac{\delta v_k}{\delta x_l} \\ & = \frac{\partial}{\partial x_l}(K_{ijkl}\frac{\delta v_{i}}{\partial x_i}\delta v_{k}) - K_{ijkl}\frac{\partial^2 \delta v_i}{\partial x_jx_l} \delta v_k \end{aligned}

Contribution 1. Energy optimization

\begin{aligned} \text{find}\quad & v(x) \\ \text{s.t.}\quad & \int_{\Omega_+} \nabla\cdot \sigma(v) dV + \int_{\Omega_+} \rho g dV = \int_{\partial \Omega_+} \lambda \cdot \hat{n}(v) dA \end{aligned}

We turn this into an optimization problem whose optimality conditions imply force balance.

\begin{aligned} \mathcal{L} & = \int_{\Omega_+} \bigg[\Psi(v) + \rho g v - \lambda(x) \phi(x + v)\bigg] dV \\ \delta\mathcal{L} & = \int_{\Omega_+}\bigg[\nabla \cdot \sigma + \rho g - \lambda(x)\cdot\hat{n}(x+v)\bigg]dV \\ & = \int_{\Omega_+}\bigg[\nabla \cdot \sigma + \rho g\bigg]dV - \int_{\partial \Omega_+} \lambda(x)\cdot \hat{n}(x+v)dA \end{aligned}

Stationarity

\begin{aligned} \text{min}_{v(x)}\quad & \int_{\Omega_+}\Psi(v)dV + \int_{\Omega_+}\rho g v dV \\ \text{s.t.}\quad & \phi_{\Omega_+}(x+v) \geq 0 \quad \forall x\in\Omega \end{aligned}

\begin{aligned} \phi(x+v) & \geq 0 \quad \forall x\\ \lambda(x) & \geq 0 \quad \forall x\\ \phi(x+v)\cdot \lambda(x) & = 0 \quad \forall x \end{aligned}

Feasibility and complementarity

Energy optimization

\begin{aligned} \text{min}_{v(x)}\quad & \int_{\Omega_+}\Psi(v)dV + \int_{\Omega_+}\rho g v dV + \int_{\Omega_+} \frac{1}{2}\rho v^2 dV\\ \text{s.t.}\quad & \phi_{\Omega_+}(x+v) \geq 0 \quad \forall x\in\Omega \end{aligned}

Generalizes quasistatic rigid-body models to the continuum mechanics setting
Contact forces (traction) is naturally handled by the non-penetration constraint
Parallel story for implicit time-stepping vs. penalty method
Solvable with Newton's method based on MPM discretization

Finally, add in a kinetic energy term to regularize the cost term

But has many drawbacks which makes me want to look for a different approach.

- Still quite costly to solve, computationally not much of a gain with second-order

- Gap in equivalence of force balance & potential minimization in presence of friction

- Does not have modeling power for plasticity.

An alternative approach

\begin{aligned} \text{min}_{v(x)}\quad & \frac{1}{2}\int_{\Omega_+}\rho v^2 dV \\ \quad & \phi_{\Omega_+}(x + v) \geq 0 \\ \quad & \textbf{vol}(\Omega_+) = \textbf{vol}(\Omega) \end{aligned}

Motivated more on position-based dynamics, with an explicit volume constraint

Intuitive physics??

Deformables in Drake

We are hoping to benchmark our method against Drake's FEM implementation.

Drake's FEM implementation:

Supports contact/deformable-rigid body coupling with SAP solver.
Only handles elastic deformation.

\Omega_+=\{x_+ | x_+= x + v(x),\; x\in\Omega\}

Deformed configuration

Exporting Deformed Configuration

We plan to export the point cloud of the dough ball in the deformed configuration and compare the positions of the points to the results obtained using our model.

--E=5000 --nu=0.4 --density=800 --time_step=1e-2

FEM Solver Failures

Experimenting with different parameters and situations that lead to solver failures:

E = 1000
nu = 0.4
density = 800
dt = 1e-2

Linear solve did not converge in Newton iterations in FemSolver. abort: Failure at multibody/contact_solvers/sap/sap_model.cc:391 in CalcDelassusDiagonalApproximation(): condition 'A_ldlt[c].eigen_linear_solver().isPositive()' failed.

E=1000 instead of 5000 and trying to squeeze the grippers closed more (closed_width = kL * 0.1 instead of kL * 0.4), it fails mid squeeze:

When time step = 1e-3 it doesn’t fail at the same place (initial gripping of the torus), but it does fail at the end (collision with ground).

E = 1000
nu = 0.4
density = 800
dt = 1e-3

FEM Solver Failures (2)

--E=5000

--nu=0.4

--density=800

--time_step=1e-2

Dropping the torus with default settings from a higher height (0.28 instead of 0.18), it would fail when it hit the ground: (you can see that the torus has penetrated the ground quite a bit)

Can be fixed by setting the time step to be 1e-3 instead of 1e-2.

--E=5000

--nu=0.4

--density=800

--time_step=1e-3

Moving Least Square Material Point Method (MLS-MPM)

"PlasticineLab: A Soft-Body Manipulation Benchmark with Differentiable Physics"

Plasticine (von Mises yield Criterion)
One-way coupling between rigid objects and soft bodies (but Compatible Particle-In-Cell (CPIC) algorithm can handle two-way coupling)
OpenAI Gym API/Environment

\Omega_+=\{x_+ | x_+= x + v(x),\; x\in\Omega\}

Deformed configuration

Key parameters & Runtimes

Key parameters:

number of particles = 10,000
grid size = 64 x 64 x 64
grid spacing = 1/n_grid = 0.0156
dt = 1e-4
E = 5000
nu = 0.2
yield stress = 50
ground friction = 100
manipulator friction = 0.9

Runtimes on 2019 Macbook Pro CPU

Taichi-Lang supports Mac GPU (Metal), but this particular implementation results in a segfault.
This video is 570 iterations of MPM
- With rendering: , about 6 minutes
- Without rendering: about 18 seconds

MLS-MPM Failure Cases

dt = 1e-3 (10x default)

dt = 2e-4 (2x default)

Courant–Friedrichs–Lewy Condition Comparison

dt = 1e-4 (default)

dt = 3e-4 (3x default)

Max Courant Number = 0.135471

Max Courant Number = 0.4 -> 2 -> 3.8 -> 18 -> 101156765360651704991744 ->inf

JMLR Ideas

How do we extend what we have right now?

In some sense, the question is always about "how to do policy gradient the right way". We've considered variation in terms of gradients, but there are other variations?

1. Where should we inject noise / take expectations to smooth out the irregular landscape? (dynamics, policy output, parameters)

Contact-Aware Control

What does stability and feedback control mean in the context of manipulation?