Mitigation of Cavitation Damage Erosion in Liquid Metal Spallation Targets Nov. 30 Dec. 1, 2005, ORNL Direct Numerical Simulation of Bubbly and Cavitating Flows and Applications to Cavitation Mitigation Roman Samulyak, Tianshi Lu In collaboration with James Glimm, Zhiliang Xu Computational Science Center Brookhaven National Lab Upton, NY 11973 Brookhaven Science Associates U.S. Department of Energy 1 Talk outline Main ideas of front tracking and the FronTier code Direct numerical simulation of cavitating and bubbly flows

Discrete vapor bubble model Dynamic bubble insertion algorithms Riemann solution for the phase boundary Validation of models Simulation of multiphase flows in the following applications: Atomization of a high speed jet Neutrino Factory/Muon Collider target Cavitation mitigation in the SNS target Conclusions and Future Plans Brookhaven Science Associates U.S. Department of Energy 2 Main ideas of front tracking Front Tracking: A hybrid of Eulerian and Lagrange methods Major components: 1) Front propagation, 2) Wave (smooth region) solution

Two separate grids to describe the solution: 1. A volume filling rectangular mesh 2. A unstructured (N-1) dimensional Lagrangian mesh to represent interface Advantages of explicit interface tracking: Real physics models for interface propagation Different physics / numerical approximations in domains separated by interfaces No interfacial diffusion Brookhaven Science Associates U.S. Department of Energy 3 The FronTier Code FronTier is a parallel 3D multiphysics code based on front tracking Physics models include Compressible fluid dynamics MHD Flow in porous media Elasto-plastic deformations Phase transition models Exact and approximate Riemann solvers, realistic EOS models

Adaptive mesh refinement Resolving interface tangling by using the grid based method Brookhaven Science Associates U.S. Department of Energy 4 Main FronTier applications Rayleigh Taylor and Richtmyer-Meshkov fluid instabilities Targets for future accelerators Tokamak refueling through the ablation of frozen deuterium pellets Brookhaven Science Associates U.S. Department of Energy Liquid jet

breakup and atomization 5 Supernova explosion Modeling of Bubbly and Cavitating Flows using the Method of Front Tracking Brookhaven Science Associates U.S. Department of Energy 6 Two models for cavitating and bubbly fluids Heterogeneous method (Discrete Bubble Model): Each individual bubble is explicitly resolved using FronTier interface tracking technique. Stiffened Polytropic EOS for liquid Polytropic EOS for

gas (vapor) Homogeneous EOS model. The mixture of liquid and vapor is treated as a pseudofluid (single-component flow); Suitable averaging is performed over a length scale of several bubbles. Small spatial scales are not resolved. Brookhaven Science Associates U.S. Department of Energy 7 Features of the discrete vapor bubble model First principles simulation. Accurate description of multiphase systems limited only by numerical errors. Resolves small spatial scales of the multiphase system.

Accurate treatment of drag, surface tension, viscous, and thermal effects. Mass transfer due to phase transition (Riemann problem for the phase boundary). Brookhaven Science Associates U.S. Department of Energy 8 Theory on Bubbly Flows Mass and Momentum Conservation u ( u 2 p) 0 t x 1 p u , 2 f c f t x t

f (1 ) g The Keller Equation 2 2 1 dR dR 1 dR R d 1 1 dR d R 3 1 R 2 1 1 pB p 2 3c dt dt c dt c dt l c dt dt 2 p g R 3 constant, p g pB R Brookhaven Science Associates U.S. Department of Energy 9

Theory on Bubbly Flows Dispersion Relation c: cf: B: : k 1 1 1 ( )2 2 2 c 1 i ( ) 2 cf B B low frequency sound speed sound speed in pure fluid resonant frequency damping coefficient 1 1

p ( ( 1 ) )( ) c , g f 2 2 2 c g cg f cf f

Steady State Shock Speed 1 1 a f U st2 c f 2 Pb Brookhaven Science Associates U.S. Department of Energy 10 B 1 3p R f Linear Wave Propagation fB R 0.06mm, 0.02%

Dispersion Relation k k ik Phase Velocity 300 2 35 theory with =0.7 simulation c in pure fluid 250 Attenuation Rate 20 log10 e k 2 theory with =0.7 simulation 30 25

(dB/cm) V (cm/ms) 1 1 3P 54.4 kHz 2R f 200 150 20 15 10 100 50 0 5 100 200

f (KHz) Brookhaven Science Associates U.S. Department of Energy 300 0 0 400 11 100 200 f (KHz) 300 400 Shock Wave Propagation

Shock speeds measured from the simulations are within 10% deviation from the steady state values. Shock Profile 2.0 The oscillation period is longer for larger bubble volume fraction. 1.8 He ( = 1.67 ) N2 ( = 1.4 ) SF6 ( = 1.09 ) 1.6 P (atm) The oscillation amplitude is smaller for gas with larger . 1.4

1.2 1.0 0 1 2 3 4 t (ms) Pa = 1.1 atm, Pb = 1.727 atm, R = 1.18 mm, = 2.5E-3 Meardured 40 cm from the interface. Brookhaven Science Associates U.S. Department of Energy 12 Shock Wave Propagation Shock profile of SF6 gas bubbles Simulation

Experiment The oscillation period is shorter than the experimental value by 28%. Brookhaven Science Associates U.S. Department of Energy [courtesy of Beylich & Glhan] 13 Dynamic Bubble Insertion Algorithm for Direct Numerical Simulation A cavitation bubble is dynamically inserted in the center of a rarefaction wave of critical strength A bubbles is dynamically destroyed when the radius becomes smaller than critical. In simulations, critical radius is determined by the numerical resolution. With AMR, it is of the same order of magnitude as physical critical radius. There is no data on the distribution of nucleation centers for mercury at the given conditions. Some estimates within the homogeneous nucleation theory: critical radius: 2S

RC PC 3 2 S W 16 S nucleation rate: J J 0 e Gb , J 0 N , Gb CR , WCR 2 m kT 3 PC 3 16 S Pc Science Associates Brookhaven U.S. Department

Energy J 0Vdt 3kTofln 1/2 Critical pressure necessary to create a 14 bubble in volume V during time dt Riemann Problem for the Phase Boundary Brookhaven Science Associates U.S. Department of Energy 15 Governing Equations and Boundary Conditions ( u ) 0 t ( u ) ( uu ) p g t ( E )

(( E p )u ) 2T t Phase Boundary Conditions (Generalized Rankine-Hugoniot Conditions): [ un ] s[ ] [ un2 p ] s[ un ] T [ Eun pu n ] s[ E ] n Brookhaven Science Associates U.S. Department of Energy 16 Interfacial Thermal Conditions 1. Equal interfacial tempareture : Tl Tv Ts The generalized Hugoniot relation : T [ ] M ev ([ H ] V [ p ]) M ev L n L : latent heat; H : Enthalphy; M ev : Mass flux V ( l v ) / 2; 1 / . 2. Two cases:

a) Contact with thermal conduction b) Phase boundary Contact with thermal conduction: Brookhaven Science Associates U.S. Department of Energy T M ev 0, ( ) 0 n 17 Phase Boundary Conditions psat (T) pv Mass flux : M ev 2RTRT : evaporation coefficient psat (T) : Clausius - Clapeyron equation pv : vapor pressure R kB ; k B is Boltzmann const.; m is molecular mass m

Brookhaven Science Associates U.S. Department of Energy 18 A deviation from Clausius-Clapeyron on vapor side is allowed. Similar to: Matsumoto etal. (94) Two Characteristic Equations t t Sl Sr Phase Boundary u c u c t

S 2 S 1 S f S 0 S 0 Sb un dp 2T c 2 d d n un dp 2T c 2 d d n Brookhaven Science Associates

U.S. Department of Energy New Position 19 S1 S2 x Phase Boundary Propagation An Iteration Algorithm. 1. Solve for mass flux and interfacial temperature by: pl pv 1 T T ( l v ) ( v v l l ) 2 M n n p (T ) pv

M sat 2RT l v 2. Solve the characteristic equations with the jump conditions: M u v , n ul , n v l p pl M 2 v v l 3. Compare the newly obtained pv and with the previous iteration to determine the convergence of the iteration. Brookhaven Science Associates U.S. Department of Energy 20 Validation Phase Boundary Solutions Brookhaven Science Associates U.S. Department of Energy

21 Applications: Liquid jet breakup and atomization Neutrino Factory / Muon Collider target Cavitation mitigation in SNS target Brookhaven Science Associates U.S. Department of Energy 22 Liquid Jet breakup and Spray Formation Breakup Regimes: 1. Rayleigh breakup 2. First wind-induced breakup 3. Second wind-induced breakup 4. Atomization DROP AND SPARY FORMATION FROM A LIQUID JET, S.P.Lin, R,D. Reitz, Annu. Rev. Fluid Mech. 1998. 30: 65-105 Brookhaven Science Associates

U.S. Department of Energy 23 Simulation setup and processes influencing atomization Inlet pressure fluctuation Cavitation in the nozzle and free surface jet Boundary rearrangement effect Brookhaven Science Associates U.S. Department of Energy 24 Simulation Results Density Plot of Jet Simulation Using Discrete Vapor Bubble Model Density Plot of Jet Simulation Using Homogenized EOS Model

Brookhaven Science Associates U.S. Department of Energy 25 Animation of Simulation Using Discrete Bubble Model Brookhaven Science Associates U.S. Department of Energy 26 Neutrino Factory / Muon Collider Target Brookhaven Science Associates U.S. Department of Energy 27 Numerical simulations of the mercury jet target Simulation of the mercury jet target interacting with a proton pulse in a magnetic field Studies of surface instabilities, jet breakup, and cavitation

MHD forces reduce both jet expansion, instabilities, and cavitation Richtmyer-Meshkov instability of the mercury target surface. Single fluid EOS (no cavitation) Brookhaven Science Associates U.S. Department of Energy 28 Cavitation in the mercury jet interacting with the proton pulse Initial density Density at 20 microseconds Initial pressure is 16 Kbar 400 microseconds Brookhaven Science Associates U.S. Department of Energy 29

MHD effects in the Mercury Target Distortion of the mercury jet by a magnetic field Brookhaven Science Associates U.S. Department of Energy Stabilizing effect of the magnetic field 30 SNS and Cavitation Mitigation P0 (r , z ) 500e r 2 0.1 z bar Courtesy of Oak Ridge National Laboratory Brookhaven Science Associates U.S. Department of Energy

31 Step 1: DNS of pressure wave propagation in the container: pure mercury Pure Mercury t Pw (t ) Pw0 e cos( Pw0 500bar 940s T 70 s Brookhaven Science Associates U.S. Department of Energy 32 2t )

T Step 2: DNS of pressure wave propagation in the container: mercury containing gas bubbles Bubbly Mercury ( R=1.0mm, =2.5% ) Pw (t ) Pw0 e t cos( Pw0 600bar 50s T 12s Brookhaven Science Associates U.S. Department of Energy 33 2t

) T Step 3: Collapse pressure of cavitation bubbles The Keller Equation 2 2 1 dR dR 1 dR R d 1 1 dR d R 3 1 R 2 1 1 pB p 2 3c dt dt c dt c dt l c dt dt 2 p g R 3 p g 0 R03 , p g pB R p(t ) P sin( 2t

0 ) T Empirical formula for P < 10Kbar and T < 1ms 1 93.0 P 1.25 p g 0 R0 3 0.50 Pc ( P, T ) Pc ( P, T , 0 0.63 ) ( ) ( ( ) ) Kbar 2 2 2 2 f cf f cf cfT Brookhaven Science Associates U.S. Department of Energy 34 Step 4: Estimation of efficiency of the cavitation damage mitigation Statistical averages of collapsing bubbles pressure peaks

Pc ( 0) E ( , R) Pc ( , R ) R0 = 1.0 m pg0 = 0.01 bar E(,R) is independent of R0 and pg0. E ~ 40 at R = 0.5 mm and 0.5% void fraction Brookhaven Science Associates U.S. Department of Energy 35 Conclusions and Future Plans Developed components enabling the direct numerical simulation of cavitating and bubbly flows Discrete vapor bubble model Dynamic bubble insertion algorithms

Riemann solution for the phase boundary Studied multiphase flows in the following applications: Atomization of a high speed jet Neutrino Factory/Muon Collider target Cavitation mitigation in the SNS target Future work: Investigate the influence of the init. bubble size on the simulation. Improve the bubble insertion algorithm - implement a conservative insertion. Studies of accelerator targets and liquid jets Brookhaven Science Associates U.S. Department of Energy 36