Rewriting the equations

APIC:

Suppose we have xp x_p and an evenly-spaced Eulerian grid xi x_{i} . Then associated with each (fixed) particle are a set of interaction weights wpi w_p^{i} that describe contribution of the value of a quantity qpnq^n_p stored at the particle position. In typical MPM literature this is given by N(x)={12∣x∣3−x2+23,0≤∣x∣<1−16∣x∣3+x2−2∣x∣+431≤∣x∣<20otherwise N(x) = \begin{cases} \frac{1}{2} |x|^3 - x^2 + \frac{2}{3}, & 0 \leq |x| < 1 \\ -\frac{1}{6} |x|^3+ x^2 - 2|x| + \frac{4}{3} & 1 \leq |x| < 2 \\ 0 & \textrm{otherwise}\end{cases} and Ni(x)=N(1h(x1−ih))N(1h(x2−jh))N(1h(x3−kh))N_i(x) = N(\frac{1}{h} \left(x_1 - ih \right)) N(\frac{1}{h} \left(x_2 - jh \right)) N(\frac{1}{h} \left(x_3 - kh \right)) which impose that points only apply to two neighbor points in each direction

We're chiefly concerned with the particle-to-grid and grid-to-particle transfers given by

min=∑pwi,pnmpDpn=∑iwi,pn(xin−xpn)(xin−xpn)⊤minvin=∑pwi,pnmp(vpn+Bpn(Dpn)−1(xin−xpn)) \begin{align*} m_i^n &= \sum_p w_{i,p}^nm_p \\ D_p^n &= \sum_i w_{i,p}^n (x_i^n - x_p^n)(x_i^n - x_p^n)^\top\\ m_i^n v_i^n &= \sum_p w_{i,p}^nm_p (v_p^n + B_p^n (D_p^n)^{-1} (x_i^n - x_p^n)) \end{align*} where BpB_p is transient and stored on particles. G2P reads as vpn+1=∑iwi,pnvin+1Bpn+1=∑iwi,pnvin+1(xin−xpn)⊤. \begin{align*} v_p^{n+1}&= \sum_i w_{i,p}^n v_i^{n+1} \\ B_p^{n+1} &= \sum_i w_{i,p}^n v_i^{n+1} (x_i^n - x_p^n)^\top. \end{align*} Together these describe a material-agnostic treatment of p2g and g2p transfers. In what follows, we make the definition Cpn=Bpn(Dpn)−1C_p^n = B_p^n (D_p^n)^{-1}.

MLS reconstruction (before quadrature)

Suppose we have samples ui=u(xi) u_i = u(x_i) . To better condition the fit, we typically choose a base point xx about which to center our linear fit, which in MPM is typically particle positions. If we have a (unconstrained) point zz at which we want to approximate, we can try to reconstruct the value based on P=[p1(z−x),…,pn(z−x)]⊤P = [p_1(z-x), \ldots, p_n(z-x)]^\top as u(z)=P(z−x)⊤c(x) u(z) = P(z-x)^\top c(x) . We choose to minimize ∑iξi(x)(P(xi−x)⊤c(x)−ui) \sum_i \xi_i(x) \left(P(x_i - x)^\top c(x) - u_i \right) where ξi\xi_i is an interaction kernel with compact support.

Solving this gives u(z)=ξi(x)P(z−x)M−1(x)P(xi−x)uiu(z) = \xi_i(x) P(z-x) M^{-1}(x)P(x_i - x)u_i with M=∑iξi(x)P(xi−x)⊤P(xi−x) M = \sum_i \xi_i(x) P(x_i-x)^\top P(x_i-x). Essentially, the second term in the sum are the weights found by MLS and the multiplication by P(z−x)⊤P(z-x)^\top follows from the equation above.

THis motivates the MLS shape functions Φi(z)=ξi(xp)P(z−xp)⊤M−1(xp)P(xi−xp) \Phi_i(z) = \xi_i(x_p)P(z - x_p)^\top M^{-1}(x_p)P(x_i - x_p)

Incompressible fixed corotated

Graphics implementations of MPM phrase constitutive relations in terms of energy functionals. The versatile one used in Stomakhin's work uses a multiplicative plasticity flow law with F=FEFP F = F_E F_P. The calculation of plastic deformation in the solid phase follows from this work where excess deformation (measured by singular values of FF) is passed into FPF_P. In the fluid phase, we set FE=(JE)1dI F_E = (J_E)^\frac{1}{d}I at the end of a time step, as deformation is effectively fully plastic.

The vanilla form of the energy functional is Ψμ(FE)=μ∥FE−RE∥F2Ψλ(det(FE))≡Ψλ(JE)=λ2(JE−1)2 \begin{align*} \Psi_\mu(F_E) &= \mu \|F_E - R_E \|_{F}^2 \\ \Psi_\lambda(\textrm{det}(F_E)) \equiv \Psi_\lambda(J_E) &= \frac{\lambda}{2}(J_E - 1)^2 \end{align*}

This gets modified to be suitable for highly incompressible materials according to Ψ^(FE)=Ψ^μ(FE)+Ψλ(JE), where Ψ^μ(FE)=Ψμ(JE−1/dFE) \hat{\Psi}(F_E) = \hat{\Psi}_{\mu}(F_E) + \Psi_\lambda(J_E), \textrm{ where } \hat{\Psi}_\mu(F_E) = \Psi_\mu(J_E^{-1/d} F_E) which eliminates any contribution of the dilational component of deformation to Ψλ(FE) \Psi_\lambda(F_E) Using concepts from hyperelasticity (see, e.g., this classic) we find σμ=1J∂Ψ^μ∂FEFE⊤ \sigma_\mu = \frac{1}{J} \pder{\hat{\Psi}_\mu}{F_E}F_E^\top . This also indicates that an analogue of the Newtonian stress-strain relationship is encoded by μ \mu . It's clear we're using the fluid-mechanical decomposition σ=σμ+σλ \sigma = \sigma_\mu + \sigma_\lambda, so we get by the chain rule σλ=1J(∂Ψλ∂JE∂JE∂FE)FE⊤=1JEJP∂Ψλ∂JEJEFE−⊤FE⊤=(1JE∂Ψλ∂JE)I≡−pI \sigma_\lambda = \frac{1}{J} \left(\pder{\Psi_\lambda}{J_E} \pder{J_E}{F_E} \right)F_E^\top = \frac{1}{J_EJ_P} \pder{\Psi_\lambda}{J_E} J_E F_E^{-\top}F_E^\top = \left(\frac{1}{J_E} \pder{\Psi_\lambda}{J_E}\right) I \equiv -pI Deviatoric stress is dealt with either implicitly or explicitly using the MLS-MPM treatment of v∗−vnΔt=1ρn∇⋅σμ+g \frac{v^* - v^n}{\Delta t} = \frac{1}{\rho_n} \nabla \cdot \sigma_\mu + g , resulting in gridpoint values of v∗v^*. In vanilla MPM, quantities needed to solve for pressure are then transfered to cell-centered grid points analogously to mass, e.g. [JE]cn=1mcn∑pwc,pnmp[JE]pn. [J_{E}]_c^n = \frac{1}{m_c^n} \sum_p w_{c, p}^n m_p [J_{E}]_p^n. This is the most important equation for determining a MLS-consistent particle-to-grid for fluids.

MLS digression

The continuum momentum equation we're solving looks like ρDvDt=∇⋅σμ+∇⋅σλ+ρg=∇⋅σμ−∇p+ρg \rho \frac{D v}{Dt} = \nabla \cdot \sigma_\mu + \nabla \cdot \sigma_\lambda + \rho g = \nabla \cdot \sigma_\mu - \nabla p + \rho g The weak form of this equation looks like 1Δt∫Ωnρ(vα∗−vαn)qα dx=∫∂ΩnTαqα ds−∫Ωnσαβ∂qα∂xβ dx\frac{1}{\Delta t} \int_{\Omega^n} \rho (v^*_\alpha - v^n_\alpha) q_{\alpha} \intd{x} = \int_{\partial \Omega^n} T_\alpha q_\alpha \intd{s} - \int_{\Omega^n} \sigma_{\alpha \beta} \pder{q_\alpha}{x_{\beta}} \intd{x}

Using the fact that ∇⋅(σq)=q⋅(∇⋅σ)+∇q:σ \nabla \cdot (\sigma q) = q \cdot (\nabla \cdot \sigma) + \nabla q : \sigma The derivation looks like ρdvdt=∇⋅σ+ρg  ⟹  ⟨ρdvdt,q⟩=⟨∇⋅σ,q⟩+⟨ρg,q⟩  ⟹  ∫Ω⟨ρdvdt,q⟩ dx=∫Ω⟨∇⋅σ,q⟩ dx+∫Ω⟨ρg,q⟩ dx  ⟹  ∫Ω⟨ρdvdt,q⟩ dx=∫Ω∇⋅(σq) dx−∫Ω∇q:σ dx+∫Ω⟨ρg,q⟩ dx  ⟹  ∫Ω⟨ρdvdt,q⟩ dx=∫∂Ω⟨σq,n⟩ ds−∫Ω∇q:σ dx+∫Ω⟨ρg,q⟩ dx \begin{align*} &\rho \der{v}{t} = \nabla \cdot \sigma + \rho g \\ \implies& \left\langle \rho \der{v}{t}, q \right\rangle = \left\langle\nabla \cdot \sigma, q \right\rangle + \left\langle\rho g, q \right\rangle\\ \implies& \int_\Omega \left\langle \rho \der{v}{t}, q \right\rangle \intd{x} = \int_\Omega \left\langle\nabla \cdot \sigma, q \right\rangle \intd{x} + \int_{\Omega} \left\langle\rho g, q \right\rangle \intd{x}\\ \implies& \int_\Omega \left\langle \rho \der{v}{t}, q \right\rangle \intd{x} = \int_\Omega \nabla \cdot (\sigma q) \intd{x} - \int_\Omega \nabla q : \sigma \intd{x} + \int_{\Omega} \left\langle\rho g, q \right\rangle \intd{x}\\ \implies& \int_\Omega \left\langle \rho \der{v}{t}, q \right\rangle \intd{x} = \int_{\partial \Omega} \left\langle \sigma q, n \right\rangle \intd{s} - \int_\Omega \nabla q : \sigma \intd{x} + \int_{\Omega} \left\langle\rho g, q \right\rangle \intd{x}\\ \end{align*} Using the MLS shape functions qα=P⊤(x−xpn)M−1(xpn)P(x−xpn)q_\alpha = P^\top(x-x_p^n) M^{-1}(x_p^n )P(x-x_p^n) and using linear polynomials for PP (APIC) and tensor-spline weights, we get ∂qα∂xβ=∂P⊤∂xβ(x−xpn)M−1(xpn)P(xi−xpn)=Mp−1Ni(xpn)(xi,β−xp,β) \begin{align*} \pder{q_\alpha}{x_\beta} &= \pder{P^\top}{x_\beta}(x-x_p^n) M^{-1}(x_p^n ) P(x_i-x_p^n) \\ &= M_p^{-1} N_i(x_p^n) (x_{i, \beta} - x_{p, \beta}) \end{align*}

Grid-space pressure correction:

The most economical grid-space equation for this (semi-incompressible) constitutive equation is JPnλnJEnpn+1Δt−δt∇⋅(1ρn∇pn+1)=JPnλnJEnpnΔt−∇⋅v∗=JPnλnJEn⋅(−1JPnλn(JEn−1))−∇⋅v∗ \begin{align*} \frac{J_P^n}{\lambda^n J_E^n} \frac{p^{n+1}}{\Delta t} - \delta t \nabla \cdot \left(\frac{1}{\rho_n} \nabla p^{n+1}\right) &= \frac{J_P^n}{\lambda^n J_E^n} \frac{p^n }{\Delta t} - \nabla \cdot v^* \\ &= \frac{J_P^n}{\lambda^n J_E^n} \cdot \left(- \frac{1}{J_P^n} \lambda^n (J_E^n - 1) \right) - \nabla \cdot v^* \end{align*} which then gives a pressure correction vn+1=v∗+Δtρn∇pn+1v^{n+1} = v^* + \frac{\Delta t}{\rho^n }\nabla p^{n+1} .

Parent post: