MPM Snow

LeemHyungJun·2023년 4월 14일
0

Computer Graphics

목록 보기
2/2
post-thumbnail

Material Point Method for Snow Simulation (SIGGRAPH 2013)


mpm........

1. Introduction

3. Paper Overview

4. Material Point Method

4-1. Full Method

1. Rasterize particle data to the grid

  • transfer mass and velocity to the grid
  • mass
    min=pmpwipnm_i^n = \sum_p m_pw_{ip}^n
  • velocity (운동량 보존을 위한 식)
    vin=pvpnmpwipn/minv_i^n = \sum_p v_p^nm_pw_{ip}^n/m_i^n

2. Compute Particle volumes and densities

  • 첫 time step 에서만 실행
  • particle의 부피 설정
  • Cell의 밀도
    mi0/h3m_i^0/h^3
  • particle의 가중치와 Volume
    ρp0=imi0wip0/h3\rho_p^0=\sum_i m_i^0w_{ip}^0/h^3
    Vp0=mp/ρp0V_p^0=m_p/\rho_p^0

3. Compute grid forces

  • Eq.6

4. Update velocities on grid

  • Section 10 using viv_i^*

5. Grid-based body collisions

  • Section 8 (Body Collision) using viv_i^*

6. Solve the linear system

  • implicit) Eq.9 for semi-implicit integration
  • explicit) vin+1=viv_i^{n+1} = v_i^* 로 계산

7. Update deformation gradient

  • Fpn+1=(I+Δtvpn+1)FpnF_p^{n+1} = (I+\Delta t\nabla v_p^{n+1})F_p^n
  • vpn+1=ivin+1(wipn)T\nabla v_p^{n+1}=\sum_i v_i^{n+1}(\nabla w_{ip}^n)^T
  • Section 7

8. Update particle velocities

  • new particle velocities
    vpn+1=(1α)vPICpn+1+αvFLIPpn+1v_p^{n+1}=(1-\alpha)v_{PICp}^{n+1}+\alpha v_{FLIPp}^{n+1}
  • PIC part
    vPICpn+1=ivin+1wipnv_{PICp}^{n+1} = \sum_i v_i^{n+1}w_{ip}^n
  • FLIP part
    vFLIPpn+1=vpn+i(vin+1vin)wipnv_{FLIPp}^{n+1}=v_p^n+\sum_i(v_i^{n+1}-v_i^n)w_{ip}^n
  • α\alpha
    α=0.95\alpha = 0.95

9. Particle-based body collisions

  • Section 8 using vpn+1v_p^{n+1}

10. Update particle positions

  • xpn+1=xpn+Δtvpn+1x_p^{n+1}=x_p^n+\Delta t v_p^{n+1}

5. Comstitutive Model

  • hyperelasticity 에 있어서 finite-strain multiplicative plasticity 를 간소화함
  • plastic yield criteria : using principal stretches
  • hyperelasticity energy density : using Lame parameters
  • Elasto-plastic energy density function
    Ψ(FE,FP)=μ(FP)FEREF2+λ(FP2(JE1)2\Psi(F_E, F_P) = \mu(F_P)||F_E- R_E||^2_F + \frac{\lambda(F_P}{2}(J_E-1)^2 - Eq(1)
    FF: Deformation Gradient (F=FEFP)F=F_EF_P)
    FEF_E : elastic part , FPF_P : plastic part
  • plastic deformation gradient
    μ(FP)=μ0eξ(1JP)\mu(F_P) = \mu_0e^{\xi(1-J_P)} and λ(FP)=λ0eξ(1JP)\lambda(F_P) = \lambda_0e^{\xi(1-J_P)} - Eq(2)
    JE=def(FE)J_E = def(F_E) and JP=det(FP)J_P = det(F_P), FE=RESEF_E = R_ES_E
    λ0,μ0\lambda_0, \mu_0 : initial Lame coefficient
    ξ\xi : dimensionless plastic hardening parameter
    θc\theta_c and θs\theta_s : compression and stretch thresholds to start plastic deformation(fracture)
    -> FEF_E의 값을 [1θc, 1+θs][1-\theta_c, \ 1+\theta_s]로 제한
  • deforming plastically
    • small deformation 영역에서는 elastic (FEF_E dependency)
    • deformation gradient 가 threshold를 넘은 경우, deforming plastically (Section.7)
    • Eq.2에 나온 material property 를 변화시킴
      • 압력(packing)에 있어서는 크게
      • stretch(farcture)에 있어서는 작게
  • Threshold
    • θc\theta_c and θs\theta_s : material이 breaking 되는 시간 결정
  • Parameters

6. Stress-based forces and linearization

  • Total elastic potential energy -> energy density Ψ\Psi
    Ω0Ψ(FE(X), FP(X))dX\int_{\Omega_0}\Psi(F_E(X),\ F_P(X))dX -Eq (3)
  • stress-based forces에 대한 MPM 공간 이산화 = Eulerian grid node material positions에 대한 에너지의 이산 근사값을 미분한 것
    -> BUT Eulerian 그리드는 변형되지 않으므로, 그리드 노드의 위치 변화는 그리드 노드의 속도에 의해 결정된다.
  • 수식
    • 변경된 그리드 노드의 위치
      xix_i : 그리드 노드 ii의 위치
      x^i=xi+Δtvi\hat{x}_i=x_i+\Delta tv_i : 현재 그리드 속도 viv_i에 대한 그리드 노드의 deformation location
    • total elastic potential
      Φ(x^)=pVp0Ψ(F^Ep(x^),FPpn)\Phi(\hat{x}) = \sum_pV_p^0\Psi(\hat{F}_{E_p}(\hat{x}),F_{P_p}^n) : total elastic potential
      (vector of all nodes x^i\hat{x}_i denote to x^\hat{x})
      VP0V_P^0 : 처음에 particle p가 차지한 volume
      FPpnF_{P_p}^n : time tnt^n에서의 particle pp의 plastic part of FF
    • F^Ep\hat{F}_{E_p}
      F^Ep\hat{F}_{E_p} : x^\hat{x}의 elastic part
      -> F^Ep(x^)=(I+i(xi^xi)(wipn)T)FEpn\hat{F}_{E_p}(\hat{x})=(I+\sum_i(\hat{x_i}-x_i)(\nabla w_{ip}^n)^T)F_{E_p}^n - Eq(4)
    • stress-based forces에 대한 MPM 공간 이산화
      fi(x^)=Φxi^=pVp0ΨFE(F^EP(x^),FPPn)(FEPn)Twipn-f_i(\hat{x})=\frac{\partial \Phi}{\partial \hat{x_i}}=\sum_pV_p^0\frac{\partial \Psi}{\partial F_E}(\hat{F}_{E_P}(\hat{x}),F_{P_P}^n)(F_{E_P}^n)^T\nabla w_{ip}^n - Eq(5)
      fi(x^)f_i(\hat{x}) : 그리드 노드 ii에서 elastic stress에 대한 힘
      -> Cauchy stress로 표현
      σp=1JpnΨFE(F^EP,FPPn)(FEPn)T\sigma_p=\frac{1}{J_p^n}\frac{\partial \Psi}{\partial F_E}(\hat{F}_{E_P},F_{P_P}^n)(F_{E_P}^n)^T 사용해서
      fi(x^)=pVpnσpwipnf_i(\hat{x})=-\sum_pV_p^n\sigma_p\nabla w_{ip}^n - Eq(6)
      Vpn=JpnVp0V_p^n=J_p^nV_p^0 : 시간 tnt^n 에 particle pp에 의해 차지된 물질의 부피
  • Grid velocity update가 목표
    • 앞에서 말한 것들은 MPM 공간 이산화와 elastic potential 간의 관계
    • Full method 3에서 말한 compute grid forces 단계에 사용
    • potential에 대한 Hessian 사용해서 elastic part의 implicit step update표현
      δfi=j2Φxi^xj^(x^)δuj=pVp0Ap(FEPn)Twipn-\delta f_i=\sum_j\frac{\partial^2\Phi}{\partial\hat{x_i}\partial\hat{x_j}}(\hat{x})\delta u_j = \sum_pV_p^0A_p(F_{E_P}^n)^T\nabla w_{ip}^n -Eq (7)
      AP=2ΨFEFE(FE(x^),FPPn):(jδuj(wjpn)TFEPn)A_P=\frac{\partial^2\Psi}{\partial F_E\partial F_E}(F_E(\hat{x}),F_{P_P}^n) : (\sum_j\delta u_j(\nabla w_{jp}^n)^TF_{E_P}^n) -Eq (8)
      -> A=C:DA=C:D 뜻은 Aij=CijklDklA_{ij}=C_{ijkl}D_{kl}

6.1 Semi-implicit update

7. Deformation gradient update

  • temporary defining
    F^PPn+1=(I+Δtvpn+1)FEPn\hat{F}_{P_P}^{n+1} = (I+\Delta t\nabla v_p^{n+1})F_{E_P}^n and F^PPn+1=FPPn\hat{F}_{P_P}^{n+1} = F_{P_P}^n
    -> 초기의 변화는 deformation gradient 의 elastic part에 영향을 미친다.
    FPn+1=(I+Δtvpn+1)FEPnFPPn=F^EPn+1F^PPn+1F_P^{n+1}=(I+\Delta t\nabla v_p^{n+1})F_{E_P}^nF_{P_P}^n=\hat{F}_{E_P}^{n+1}\hat{F}_{P_P}^{n+1}
  • Next Step
    • F^EPn+1\hat{F}_{E_P}^{n+1}중에서 critical deformation threshold를 넘는 것을 F^PPn+1\hat{F}_{P_P}^{n+1}에 넣기
    • SVD 써서 F^EPn+1=UPΣ^pVpT\hat{F}_{E_P}^{n+1}=U_P\hat{\Sigma}_pV_p^T
      Σp=clamp(Σp^,[1θc,1+θs])\Sigma_p = clamp(\hat{\Sigma_p}, [1-\theta_c, 1+\theta_s])
  • Final
    • elastic and plastic component of deformation gradient
      FEPn+1=UpΣpVpTF_{E_P}^{n+1}=U_p\Sigma_pV_p^T and FPPn+1=VpΣp1UpTFpn+1F_{P_P}^{n+1}=V_p\Sigma_p^{-1}U_p^TF_p^{n+1}
      따라서 Fpn+1=FEpn+1FPpn+1F_p^{n+1}=F_{E_p}^{n+1}F_{P_p}^{n+1}

8. Body Collisions

  • collision check twice in each time step
    • 1st : grid velocity에 힘이 가해진 직후의 grid velocity viv_i^*
    • 2nd : particle velocity vpn+1v_p^{n+1}의 위치를 업데이트 하기 전 보간에 의한 소실을 막기 위해서
  • in case of collision
    n=ϕn=\nabla \phi : local normal
    vcov_{co} : object velocity
    • vv : particle/grid velocity
      -> transform into reference frame of the collision object
      -> vrel=vvcov_{rel}=v-v_{co}
      if) bodied separate (vn=vreln0)v_n = v_{rel}\cdot n\ge 0) -> no coliision
      let) vt=vrelnvnv_t=v_{rel}-nv_n : tangential portion of the relative velocity
      -> if) sticking impluse required (vtμvn)||v_t|| \le -\mu v_n) -> vrel=0v'_{rel}=0
      -> if) dynamic friction vrel=vt+μvnvt/vtv'_{rel}=v_t+\mu v_nv_t/||v_t|| (μ\mu : coefficient of friction)
  • Finally
    v=vrel+vcov'=v'_{rel}+v_{co}
    -> relative velocity back into world coordinates
  • Types of Collision objects
    • Rigid Case
    • Deforming Case
  • Finally
    • sticky collision
      • vrel=0v'_{rel}=0

0개의 댓글