Stability Analysis of Material Point Method Dr. Martin Berzins, Dr. - PowerPoint PPT Presentation

Stability Analysis of Material Point Method Dr. Martin Berzins, Dr. Mike Kirby, Chris Gritton Scientific Computing and Imaging Institute, University of Utah March 14, 2013

Compressible Flow in One Dimension ∂ρ ∂ t + u ∂ρ ∂ x + ρ∂ u ∂ x = 0 � ∂ u � ∂ t + u ∂ u + ∂ p ρ ∂ x = 0 ∂ x p = a 2 ρ

Compressible Flow - Brackbill’s Formulation ∂ρ 1 ∂ρ 1 ∂ u 1 ∂ t + U 0 ∂ x + ρ 0 ∂ x = 0 � ∂ u 1 � ∂ u 1 + ∂ p 1 ρ 0 ∂ t + U 0 ∂ x = 0 ∂ x Substituting in p 1 = a 2 ρ 1 and rewriting the above equations in terms of their Lagrangian and Eularian parts. Note the Lagrangian particles move at constant velocity U 0 . D ρ ∂ u i = − a 2 ∂ρ i Du = − ρ 0 ∂ x , ρ 0 Dt Dt ∂ x ���� � �� � Lagrangian Lagrangian J.U. Brackbill, The Ringing Instability in Particle-in-Cell Calculations of Low-Speed Flow, Journal of Computational Physics, 75, 469-492, 1988

The Algorithms Formulation 1 Formulation 2 i = � np i = � np 1. u t p =1 u t 1. ρ t p =1 ρ t p p i = � np 2. ρ t p =1 ρ t 2. u t +1 p + c dt = u t h ( ρ t i +1 − ρ t i ) p p = � np 3. u t +1 = u t p + c dt h ( ρ t i +1 − ρ t 3. u t +1 p =1 u t +1 i ) p p i 4. ρ t +1 = ρ t p + c dt h ( u t i +1 − u t h ( u t +1 i +1 − u t +1 i ) 4. ρ t +1 = ρ t p + c dt ) p p i 5. x t +1 = x t 5. x t +1 p + vdt = x t p + vdt p p

Formulation 1 - Analysis What happens at the nodes? np � u t +1 S t +1 u t +1 = p i ip p =1

Formulation 1 - Analysis Substituting in u t +1 from the algorithm we get, p np � u t +1 S t +1 u t +1 = i ip p p =1 np p + c dt � S t +1 [ u t h ( ρ t +1 i +1 − ρ t = i )] ip p =1 np np p + c dt � � S t +1 u t S t +1 ( ρ t i +1 − ρ t = i ) ip ip h p =1 p =1 np np p + c dt � � = u t ( S t +1 − S t ip ) u t S t +1 ( ρ t i +1 − ρ t i + i ) . ip ip h p =1 p =1

Formulation 1 - Analysis If we apply the same steps to ρ t +1 and substitute we get the i following set of equations Let the Courant term be k = c dt h . np np � � u t +1 = u t ( S t +1 − S t ip ) u t S t +1 ( ρ t i +1 − ρ t i + p + k i ) i ip ip p =1 p =1 np np � � ρ t +1 ( S t +1 S t +1 = ρ t − S t ip ) ρ t ( u t i +1 − u t i + p + k i ) . i ip ip p =1 p =1

Formulation 1 - Analysis What can we say about the term np � S t +1 ( ρ t i +1 − ρ t i )? ip p =1

Formulation 1 - Analysis Remember that p + c dt u t +1 = u t h ( ρ t i +1 − ρ t i ) . p Therefore np � S t +1 ( ρ t i +1 − ρ t i ) = C 1 ( ρ t i − ρ t i − 1 ) + C 2 ( ρ t i +1 − ρ t i ) ip p =1 C 1 + C 2 = 1

Formulation 1 - Analysis If the particle velocity is zero then, C 1 = C 2 = . 5 and � np p =1 ( S t +1 − S t ip ) = 0, and we get the following. ip u t +1 = u t i + k ( ρ t i +1 − ρ t i − 1 ) i ρ t +1 = ρ t i + k ( u t i +1 − u t i − 1 ) . i Using Von Neumann analysis we can gain some insight into the stability of the underlying schemes. Let ξ and η represent the errors at the nodes for u and ρ respectively. ξ = aG t e i β j η = bG t e i β j

Formulation 1 - Analysis Expression of error growth at the nodes. aG t +1 e i β j = aG t e i β j + k ( bG t e i β j +1 − bG t e i β j − 1 ) bG t +1 e i β j = bG t e i β j + k ( aG t e i β j +1 − aG t e i β j − 1 ) After rearranging of terms. a (1 − G ) + kb ( e i β − e − i β ) = 0 b (1 − G ) + ka ( e i β − e − i β ) = 0

Formulation 1 - Analysis The system of equations. � (1 − G ) � � a � � 0 � k (2 i sin β ) = k (2 i sin β ) (1 − G ) 0 b We can use the fact that the determinant of the system is zero and set γ = k (2 i sin β ) we get the following. (1 − G ) 2 − γ 2 = 0 G 2 − 2 G + 1 − γ 2 = 0 Solving for G we get, � G = 2 ± 4 − 4(1 − γ 2 ) 2 = 1 ± γ. This is an Unstable Scheme

Formulation 2 - Analysis np np � � u t +1 = u t ( S t +1 − S t ip ) u t S t +1 ( ρ t i +1 − ρ t i + p + k i ) i ip ip p =1 p =1 np np � � ρ t +1 = ρ t ( S t +1 − S t ip ) ρ t S t +1 ( u t +1 i +1 − u t +1 i + p + k ) . i ip ip i p =1 p =1 Note the updated time step

Formulation 2 - Analysis Again as in the previous formulation set particle velocity to zeros. u t +1 = u t i + k ( ρ t i +1 − ρ t i − 1 ) i ρ t +1 = ρ t i + k ( u t +1 i +1 − u t +1 i − 1 ) . i

Formulation 2 - Analysis Using Von Neumann analysis, aG t +1 e i β j = aG t e i β j + k ( bG t e i β j +1 − bG t e i β j − 1 ) bG t +1 e i β j = bG t e i β j + k ( aG t +1 e i β j +1 − aG t +1 e i β j − 1 ) After rearranging of terms. a (1 − G ) + kb ( e i β − e − i β ) = 0 b (1 − G ) + kaG ( e i β − e − i β ) = 0

Formulation 2 - Analysis The system of equations. � � � a � � 0 � (1 − G ) k (2 i sin β ) = kG (2 i sin β ) (1 − G ) 0 b Solve the determinant of the system and set γ = k (2 i sin β )). (1 − G ) 2 − G γ 2 = 0 G 2 − G ( γ 2 + 2) + 1 = 0 Solving for G we get, � G = ( γ 2 + 2) ± ( γ 2 + 2) 2 − 4 2 � = 1 + γ 2 ± γ 4 + 2 γ 2 2

Formulation 2 - Analysis Given γ = k (2 i sin β ) what do we learn form G ? � G = 1 + γ 2 ± γ 4 + 2 γ 2 2 � 16 k 2 sin 4 β − 8 k 2 sin 2 β = 1 − 2 k 2 sin 2 β ± 2 � = 1 − 2 k 2 sin 2 β ± k sin β 4 k 2 sin 2 β − 2 1 If k ≤ √ 2 � 1 − 2 k 2 sin 2 β | G | =

Example 1

1d Linear Elastic Model Starting with the Cauchy Momentum Equation ρ Dv Dt = ∇ · σ + b Rewriting for the 1d form and substituting for σ and neglecting the body force. � � ρ Dv Dt = ∂ E ∂ u ∂ x ∂ x � Given ρ and E as constants we get, c = E /ρ , Dt = c 2 ∂ 2 u Dv ∂ x 2 .

1d Linear Elastic Model The Coupled Set of Equations forming the Wave Equation  Du  Dt = v   ⇒ ∂ 2 u ∂ t 2 = c 2 ∂ 2 u  = ∂ x 2 Dt = c 2 ∂ 2 u Dv     ∂ x 2

The Algorithm i = � np 1. u t p =1 S ip u t p i = � np 2. v t p =1 S ip v t p i = c 2 3. a t h 2 ( u t i +1 − 2 u t i + u t i − 1 ) 4. v t +1 = v t i + dt ∗ a t i i p + dt � np 5. v t +1 = v t p =1 S ip a t p i p + dt � np 6. u t +1 = u t p =1 S ip v t +1 p i

Analysis As we did in the previous formulations lets see what happens at the nodes for u i . np � u t +1 S t +1 u t +1 = . p i ip p =1

Analysis Substituting in for u t +1 from the above algorithm we get the p following. np � u t +1 S t +1 u t +1 = i ip p p =1 np np � � S t +1 ( u t S t ip v t +1 = p + dt ) ip i p =1 p =1 np np np � � � S t +1 S t +1 ip v t +1 u t S t = p + dt ip ip i p =1 p =1 p =1 np np np � � � = u t ( S t +1 − S t ip ) u t S t +1 S t ip v t +1 i + p + dt . ip ip i p =1 p =1 p =1

Analysis If the following condition holds np � S ip = 1 , p =1 then we get the following np np np � � � u t +1 = u t ( S t +1 − S t ip ) u t S t +1 S t ip v t +1 i + p + dt i ip ip i p =1 p =1 p =1 np � ( S t +1 p + v t +1 = u t − S t ip ) u t i + dt . ip i p =1

Analysis Now lets look at v i . np � v t +1 S t +1 v t +1 = p i ip p =1 np np � � S t +1 ( v t S t ip a t = p + dt i ) ip p =1 p =1 np np np � � � = v t ( S t +1 − S t ip ) v t S t +1 S t ip a t i + p + dt i ip ip p =1 p =1 p =1 np � ( S t +1 = v t − S t ip ) v t p + a t i + i dt . ip p =1

Analysis Now substitute in for a i np � v t +1 = v t ( S t +1 − S t ip ) v t p + a t i + i dt i ip p =1 np p + c 2 dt � = v t ( S t +1 − S t ip ) v t h 2 ( u t i +1 − 2 u t i + u t i + i − 1 ) ip p =1

Analysis We can now combine the two equations, u t +1 and v t +1 . i i np � u t +1 = u t ( S t +1 − S t ip ) u t p + v t +1 i + dt i ip i p =1 np np � � ( S t +1 ( S t +1 i dt 2 = u t − S t ip ) u t p + v t − S t ip ) v t p + a t i + i dt + ip ip p =1 p =1 Knowing np � u t i = u t − 1 ( S t ip − S t − 1 ) u t − 1 + v t + i dt p i ip p =1 we can arrive at v t i dt np � v t i dt = u t i − u t − 1 ( S t ip − S t − 1 ) u t − 1 − p i ip p =1

Analysis After substituting for v t i and rearranging some terms we get the following finite difference scheme at the nodes, u t +1 i + u t − 1 = c 2 u t i +1 − 2 u t i + u t − 2 u t i − 1 i i + S , dt 2 h 2 where np 1 � ( S t +1 ip − S t − 1 ) u t − 1 − S t ip )( u t p + v t p dt ) − ( S t S = . p dt 2 ip ip p =1

Analysis Von Neumann analysis gives us the following. G 2 − 2 G γ + 1 = 0 , γ = 1 − 2 k 2 sin 2 ( β/ 2) � γ 2 − 1 G = γ ± � 1 − γ 2 is real for all β and then we get, If k ≤ 1 then | γ | ≤ 1 and � 1 − γ 2 G = γ ± i � γ 2 + (1 − γ 2 ) | G | = = 1 . John Noye, Finite Difference Techniques for Partial Differential Equations, Computational Techniques for Differential Equations, Elsevier Science Publishers B.V., p.285, 1984

Example 2

Conclusions and Further Explorations 1. The underlying formulation matters when it comes to stability. 2. How does the source term contribute to stability? 3. How are ringing instability and formulation instability connected?