What’s Langrangian view and Eulerian view? How to turn a mass spring system into linear equations system, and how to solve it? What is the most advanced research on this topic?
Table of Contents
Langrangian View and Eulerian View
拉格朗日 随波逐流
Langrangian View: sensors that move passively with the simulated material
“What are my position and velocity?”
Mostly particle physics, position based. Every particle has position and velocity.
欧拉视角 岿然不动
Eulerian View: sensors that never move
“What is the material velocity passing by?”
In Eulerian View, the field like flow velocity is represented as a function of position and time.
Related link: Lagrangian and Eulerian specification of the flow field
Mass Spring System: Easiest Analysis
Mass spring system
\[\begin{align*} \vec{f_{ij}}&=-k(||\vec{x}_i-\vec{x}_j||_2 - l_{ij})(\widehat{\vec{x}_i-\vec{x}_j})\\ \vec{f}_i&=\sum^{j\neq i}_j\vec{f}_{ij}\\ \frac{\partial\vec{v}_i}{\partial t}&=\frac{1}{m_i}\vec{f_i}\\ \frac{\vec{x}_i}{\partial t}&=\vec{v}_i \end{align*}\]$k$:spring stiffness; $l_{ij}$: spring rest length between particle i and particle j;
To calculate, we have two ways to do it.
Explicit:
- Future depends only on past
- Easy to implement
-
Easy to explode:
\[\Delta t\leq c\sqrt{\frac{m}{k}}(c\sim 1)\] - Bad for stiff materials(k is too large)
The updating frequency should allow the system to get close to steady situation. Imagine a mass spring that the stiffness is too large, so the next frame the whole system already move past the steady point and get much far away in opposite direction, then the whole system will “explode”.
Implicit:
- Future depends on both future and past
- Need to solve a system of (linear) equations
- In general harder to implement
- Each step is more expensive but time steps are larger
- Numerical damping and locking.
Solving Mass Spring System
To solve the whole system, there are three approaches:
-
Forward Euler(explicit)
\[\begin{align*} \vec{v}_{t+1}&=\vec{v}_t + \Delta t\frac{\vec{f}_t}{m}\\ \vec{x}_{t+1}&=\vec{x}_t + \Delta t\vec{v}_t \end{align*}\] -
Semi-implicit Euler(aka. symplectic Euler, explicit, commonly used)
\[\begin{align*} \vec{v}_{t+1}&=\vec{v}_t + \Delta t\frac{\vec{f}_t}{m}\\ \vec{x}_{t+1}&=\vec{x}_t + \Delta t\vec{v}_{t+1} \end{align*}\] -
Backward Euler(often with Newton’s method, implicit)
$v_{t+1}$ relies to force $f(x_{t+1})$, we’ll use $M^{-1}$ as a abbreviation from position to force.
\[\begin{align*} \vec{x}_{t+1}&=\vec{x}_t + \Delta t\vec{v_{t+1}}\\ \vec{v}_{t+1}&=\vec{v}_t + \Delta t M^{-1}\vec{f}(\vec{x}_{t+1}) \end{align*}\]This is a set of equations. To eliminate $\vec{x}_{t+1}$:
\[\vec{v}_{t+1} = \vec{v}_t + \Delta t M^{-1}f(\vec{x}_t + \Delta t\vec{v}_{t+1})\]$M^{-1}$ is ordinarily not linear, so the equation is hard to solve. To linearize(one step of Newton’s method, which is a Taylor Expansion):
\[\vec{v}_{t+1}=\vec{v}_t + \Delta t M^{-1}[\vec{f}{x_t} + \frac{\partial f}{\partial x}(x_t)\Delta t \vec{v}_{t+1}]\]Clean it up, we finally get a linear equation:
\[[I-\Delta t^2 M^{-1}\frac{\partial f}{\partial x}(\vec{x}_t)] \vec{v}_{t+1} = \vec{v}_t + \Delta t M^{-1}f(\vec{x}_t)\]
Solving Linear System: Jacobi Iterations
How to solve the linear system?
- Jacobi/Gauss-Seidel iterations (easy to implement)
- Conjugate gradients (later in this course)
We’ll start with the Jacobi iterations. First, we’ll simplify the equation to a linear form:
\[\begin{align*} A&=[I-\Delta t^2 M^{-1}\frac{\partial f}{\partial x}(x_t)]\\ b&=v_t + \Delta tM^{-1}f(x_t) \end{align*}\]Now we have a general linear equation $Av_{t+1}=b$. In Jacobi iterations, A can be decomposed into the sum of a diagonal component, a lower triangular part and an upper triangular part.
\[A = D + L + U\] \[D=\begin{bmatrix} a_{11} & 0 & \cdots & 0 \\ 0 & a_{22} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & a_{nn} \end{bmatrix}\]Then we can build a ‘fake equalization’:
\[D x^{(k+1)} + (L+U)x^{(k)} = b\]Find solution in an iteration process:
\[x^{(k+1)}=D^{-1}(b-(L+U)x^{(k)})\]The minimum amount of storage is two vectors of size n. Jacobi iteration is useful because sometimes A is a sparse matrix, and $A^{-1}$ is hard to calculate, also may not be a sparse matrix anymore which requires much more storage. The process of Jacobi iteration use the original form of A without generate any other matrices.
Note: Jacobi iteration does not always converge. We can easily prove that
\[v^{(k+1)}-v=D^{-1}(L+U)(v^{(k)}-v)\]The condition of convergence is $\rho (D^{-1}(L+U))<1$. Since we are discussing a sparse matrix, this condition seems apparently satisfied. For more information, check here.
Advanced Reading
- Smoothed-particle hydrodynamics(SPH)
- Weakly Cmpressible SPH(WCSPH) to simulate water
- SPH in computer graphics(Presentation)
- SPH techniques for the physics based simulation of fluid and solids(Github Page)
- Courant-Friedrichs-Lewy condition(CFL)
- Smoothed Particle Hydrodynamics (SPH) implemented with C++ and CUDA(Github Page)
- Interlinked SPH Pressure Solvers for Strong Fluid-Rigid Coupling(C. Gissler et al.)
- 耦合：两种不同的物理模拟系统能够双向交互
- From particles to Surfaces: Marching Cube, VDB
Other partical-based simulation methods:
- Discrete element method(Particle based simulation for granular materials, 2005)
- Moving particle semi-implicit (Moving-particle semi-implicit method for fragmentation of incompressible fluid, 1996)
- Power Particles: An incompressible fluid solver based on power diagrams(Particles: an incompressible fluid solver based on power diagrams, 2015)
- A peridynamic perspective on spring-mass fracture(J. A. Levine et al.(2014))