Viscous Compaction Waves Joseph M. Powers and Michael T. Cochran - PowerPoint PPT Presentation

Viscous Compaction Waves Joseph M. Powers and Michael T. Cochran University of Notre Dame, Notre Dame, Indiana SIAM Annual Meeting New Orleans, Louisiana 12 July 2005

Compaction Wave Schematic 100 µ m D ~ 400 m/s u ~ 100 m/s p particles φ ~ 0.73 φ ~ 0.98 s s

Introduction • Heterogeneous energetic solids composed of 100 µm crystals in plastic binder. • Engineering length scales on the order of many cm . • Macrobehavior (ignition, etc.) strongly linked to microstructure. • Continuum mixture models with non-traditional constitutive theories needed to capture grain scale physics.

Review • Gokhale and Krier, Prog. Energy Combust. Sci. , 1982. • Baer and Nunziato, Int. J. Multiphase Flow , 1986. • Powers, Stewart, Krier, Combust. Flame , 1990. • Saurel and Abgrall, J. Comput. Phys. , 1999. • Bdzil, et al. Phys. Fluids , 1999, 2001. • Gonthier and Powers, J. Comput. Phys. , 2000. • Powers, Phys. Fluids , 2004.

Issues with Continuum Mixture Theories • Well-posedness not always straightforward. • Second law complicated. • Shock jumps not clearly defined for non-conservative equations. • Consequent numerical difficulties.

Inviscid Theory of Bdzil, et al. • First theory to unambiguously satisfy the second law. • Hyperbolic and well-posed for initial value problems. • Fundamentally non-conservative. • Some regularization needed for discontinuities. • No viscous cutoff mechanism for multidimensional instabilities. • Grid-dependent numerical viscosity problematic.

Viscous Extension ∂ ∂t ( ρ s φ s ) + ∇ · ( ρ s φ s u s ) = C , ∂ ∂t ( ρ g φ g ) + ∇ · ( ρ g φ g u g ) = −C , ∂ � � ρ s φ s u s u T ∂t ( ρ s φ s u s ) + ∇ · s + φ s ( p s I − τ s ) = M , � � ∂ ρ g φ g u g u T ∂t ( ρ g φ g u g ) + ∇ · g + φ g ( p g I − τ g ) = −M ,

Viscous Extension (cont.) � � �� ∂ e s + 1 ρ s φ s 2 u s · u s ∂t � � � � e s + 1 + ∇ · ρ s φ s u s + φ s u s · ( p s I − τ s ) + φ s q s = E , 2 u s · u s � � �� ∂ e g + 1 ρ g φ g 2 u g · u g ∂t � � � � e g + 1 + ∇ · ρ g φ g u g + φ g u g · ( p g I − τ g ) + φ g q g = −E , 2 u g · u g ∂ρ s ∂t + ∇ · ( ρ s u s ) = − ρ s F φ s , ∂ ∂t ( ρ s φ s η s + ρ g φ g η g ) � φ s q s � + φ g q g + ∇ · ( ρ s φ s u s η s + ρ g φ g u g η g ) ≥ −∇ · . T s T g

Constitutive Equations φ g + φ s = 1 , ψ s = ˆ ψ s ( ρ s , T s ) + B ( φ s ) , ψ g = ψ g ( ρ g , T g ) , � � ∂ψ s ∂ψ g � � p s = ρ 2 p g = ρ 2 , , � � s g ∂ρ s ∂ρ g � � T s ,φ s T g � � η s = − ∂ψ s η g = − ∂ψ g � � , , � � ∂T s ∂T g � � ρ s ,φ s ρ g � β s = ρ s φ s ∂ψ s � , � ∂φ s � ρ s ,T s e s = ψ s + T s η s , e g = ψ g + T g η g ,

Constitutive Equations (cont.) � ( ∇ u s ) T + ∇ u s � − 1 τ s = 2 µ s 3( ∇ · u s ) I , 2 � ( ∇ u g ) T + ∇ u g � − 1 τ g = 2 µ s 3( ∇ · u g ) I , 2 q s = − k s ∇ T s , q g = − k g ∇ T g , C = C ( ρ s , ρ g , T s , T g , φ s ) , M = p g ∇ φ s − δ ( u s − u g ) + 1 2( u s + u g ) C , � � e s − u s · u s E = H ( T g − T s ) − p g F + u s · M + C , 2 F = φ s φ g µ c ( p s − β s − p g ) .

Equations of State Modified Tait equation for solid (correction courtesy D. W. Schwendeman) � T s � ρ s � � �� + 1 ρ s 0 ψ s ( ρ s , T s , φ s ) = c vs T s 1 − ln + ( γ s − 1) ln ρ s ε s + q T s 0 ρ s 0 γ s   � 2 − φ s 0 � (1 − φ s ) 1 − φs +( p s 0 − p g 0 ) (2 − φ s 0 ) 2 2 − φs ln   � � 1 − φs 0 2 − φ s 1 ρ s 0 φ s 0 ln (1 − φ s 0 ) 2 − φs 0 1 − φ s 0 Virial equation for gas � T g � ρ g � � � � �� ψ g ( ρ g , T g ) = c vg T g 1 − ln + ( γ g − 1) ln + b g ( ρ g − ρ g 0 ) T g 0 ρ g 0

Viscous Dissipation Function   ∇ u s + ( ∇ u s ) T 1   Φ s = 2 µ s − 3( ∇ · u s ) : I   2   � �� � � �� � strain rate mean strain rate � �� � deviatoric strain rate   ∇ u s + ( ∇ u s ) T 1   − 3( ∇ · u s ) . I   2   � �� � � �� � strain rate mean strain rate � �� � deviatoric strain rate • similar expression for Φ g .

Dissipation: Clausius-Duhem Equation � β s � ρ s T s + e s − e g − p g (1 /ρ g − 1 /ρ s ) I ≡ ( −C ) + η g − η s T g + δ ( u s − u g ) · ( u s − u g ) T g + H ( T g − T s ) 2 T g T s ( p s − β s − p g ) 2 + φ s φ g µ c T s + φ s Φ s + φ g Φ g T s T g + k s φ s ∇ T s · ∇ T s + k g φ g ∇ T g · ∇ T g ≥ 0 . T 2 T 2 s g

Characteristics • Three real characteristics u s , u s , u g , • Three associated eigenvectors, • Not enough eigenvectors for eleven equations: parabolic, • Eight additional conditions from boundary conditions on T s , T g , u s , u g .

Numerical Method • One-Dimensional: Fortran 90 code – Second order central spatial discretization – High order implicit integration in time with DLSODE • Two-Dimensional: FEMLAB software tool – Finite element method for the form ∂ q ∂t + ∇ · f ( q ) = s ( q ) . – Unstructured mesh

1D Verification: Shock Tube T (K) inviscid L (K) 310 T g 1 analytical 1 slope=0.75 viscous 305 T s numerical 0.1 300 0.01 slope=1.95 295 A2 A1 290 0.001 ∆ x (m) 0.5 x (m) 0.1 0.2 0.3 0.4 0.0001 0.001 0.01 0.1

1D Verification: Piston-Driven Shock T (K) 340 solid 330 gas 320 310 B1 300 x (m) 0.1 0.2 0.3 0.4 0.5 T (K) T (K) g s viscous shock in solid viscous shock in gas 340 325 steady solution steady solution 330 320 time-dependent time-dependent 315 solution 320 solution 310 B2 310 B3 305 300 0.39 x (m) 0.28 x (m) 0.36 0.37 0.38 0.18 0.20 0.22 0.24 0.26

1D Subsonic Piston-Driven Compaction φ s p (MPa) 1 50 p β s s 0.8 10 0.6 5 0.4 p 1 g E1 0.2 E2 0.5 x (m) x (m) 0 0.1 0.2 0.3 0.4 0.5 0.1 0.2 0.3 0.4 0.5 T (K) u (m/s) 304 100 80 303 T , T u , u g s s g 60 302 40 301 E3 20 0.5 x (m) x (m) 0.1 0.2 0.3 0.4 0.1 0.2 0.3 0.4 0.5

1D Dissipation: Subsonic Case 3 I (MW/m /K) 200 Total Compaction 150 100 E 50 Solid Momentum Diffusion x (m) 0.1 0.2 0.3 0.4 0.5

FEMLAB vs. F90 Verification: 1D Shock Tube 1.4 x 10 -3 310 305 1.0 r o r r E T [K] e 300 v i g t a l 0.5 e R 295 0 290 0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5 x [m] x [m]

Narrow 2D Shock Tube vs. 1D Shock Tube 300 K 292 K 310 K 300 K 2 D 0 0.1 0.2 0.3 0.4 0.5 x (m) 310 1 D t = 50 µ s T (K) 300 g 290 0.0 0.1 0.2 0.3 0.4 0.5 x (m)

Small Energy Pulse: 2D Response . t = 18 μ s small temperature perturbation . radial pressure wave at ~ 2000 m/s . t = 0 μ s reflection at wall Solid Pressure

Large Energy Pulse: 2D Response Max = 0.867 φ s (x,y) 0.040 y (m) 0.000 -0.015 0 0.035 0.05 Min=0.730 x (m)

Conclusions • Diffusion enables use of simple numerical techniques. • Diffusion suppresses short wavelength instabilities, e.g. Kelvin-Helmholtz. • Diffusion suppresses subgranular length scales. • Compaction dominates the dissipation. • Rigorous subscale physical justification for diffusion models presently lacking. • Such justification necessary for a validated model.