2D simulations of the dynamics and fracture of metals in the energy absorption zone of the high-current electron beam

The metal behavior under the action of high-current electron irradiation is an example of a complex problem with a large number of involved physical processes. In this paper, we develop a continuum model of dynamic fracture of a metal melt in the energy absorption zone of a high-current electron beam and present a two-dimensional version of this model. The 2D model is necessary for the description of the action of a narrow-focused electron beam on matter. Fracture of copper and aluminum in the molten state in the energy absorption zone of the electron beam is numerically simulated, and the spatial distribution of the droplet size is investigated. Action of a high-current electron beam with the electron energy of about 1 MeV leads to the formation and subsequent breakup of the melt to droplets with diameters from 0.5 μm to several tens of micrometers. The finest droplets (0.5-5 μm) are generated in the center of the energy absorption zone, while the droplet sizes increase up to several tens of micrometer at the edge of the energy absorption zone, where matter is colder. Keeping the melt in a metastable (expanded and overheated) state provides a tension wave with an amplitude of about 3-4 GPa following the shock wave and propagating deep into the target.


Introduction
High-current electron irradiation [1][2][3][4][5][6] induces a large number of different physical processes in irradiated matter. An intensive energy release brings matter into an extreme state with high values of pressure and temperature. Leaving the extreme state is accompanied by a number of dynamic processes including propagation of compression and tension waves [1,3,6], fracture in the solid state [6] due to generation and growth of voids, melting and fracture in the molten state [7][8][9][10] due to cavitation, evolution of the two-phase medium consisting of liquid and vapor, and melt fragmentation to liquid droplets [11]. In the present study, we consider the evolution of the molten state of metal, in particular, the processes of cavitation, evaporation and fracture in the molten state.
Dynamic fracture of liquids is most often investigated using theoretical and numerical methods. The number of experimental works on melt fracture [12,13] is substantially restricted by the high complexity of parameters diagnostics for matter in the extreme state. Molecular dynamics is widely and fruitfully used [14][15][16], especially for the investigation of processes in matter under the action of ultra-short intensive pulses of laser irradiation. In the case of highcurrent electron irradiation, the length and time scales are much larger than those accessible by molecular dynamics. Therefore, a continuum model should be developed and used for the description of high-current electron irradiation. In this paper, we propose and use a 2D variant of a previously developed 1D model of fracture of a metal melt.

Mathematical model
We consider a 2D cylindrical geometry of the problem, which allows us to investigate normal incidence of a focused high-current electron beam. In this section, we give the main equations.
The evolution of metal under irradiation can be divided into the following three stages: S1. Homogeneous condensed liquid or solid metal. S2. Melt (liquid metal) with vapor bubbles. S3. Vapor with droplets of liquid phase.
Transition from S1 to S2 is connected with thermal nucleation of supercritical bubbles (cavities) in the metastable (expanded or overheated) liquid phase (cavitation) [14,15,17]. At stage S2, the substance is a bubbly liquid [18]. Transition from S2 to S3 includes the fragmentation of the liquid phase to separate droplets, which takes place at a high enough volume fraction of the vapor phase. The liquid droplets tend to become spherical due to the action of surface tension. At stage S3, the vapor fills a singly connected area.
The dynamics of matter is described by two-phase continuum mechanics equations [19]. An approximation of interpenetrating continuums is used: both phases are described by continuous fields of parameters, such as concentrations and radii of bubbles or droplets, pressures, temperatures, densities and volume fractions of liquid and vapor phases. The solid or liquid metal (melt) is a singly connected carrying agent at stages S1 and S2, while the vapor is the carrying agent at stage S3. The vapor bubbles form a dispersed phase at stage S2, while the liquid droplets are the dispersed phase at stage S3. A one-velocity approximation is used as well: velocities of the carrying agent and the dispersed phase are assumed to be equal.

Conservation laws
In the framework of the one-velocity approximation, the equations of mechanics for mixture as a whole are similar to the case of a one-phase medium. The Lagrangian frame of reference is used. Velocity is determined by stresses in the carrying agent [19], and the equation of motion takes the following form: where r and z are cylindrical coordinates; v r and v z are the components of the substance velocity; ρ is the average mass density of the mixture; P c is the pressure in the carrying agent; S rr , S zz and S rz are the components of the tensor of stress deviators in the carrying agent. These stress deviators are non-zero only in the case of a solid carrying agent, and they are calculated using the dislocation plasticity model [20,21]. The average density of mixture is connected with the true densities of the carrying agent ρ c and the dispersed phase ρ d by the expression ρ = ρ c α c + ρ d α d , where α c and α d are the volume fractions of the corresponding phases. The subscript "c" denotes the values characterizing the carrying agent, and the subscript "d" relates to the dispersed phase. The equation of continuity for mixture takes the following form: where U = (U c α c ρ c + U d α d ρ d ) /ρ is the specific average internal energy of the mixture, U c and U d are the specific internal energies of the carrying agent and dispersed phases, respectively. In Eq. (4): D is the energy release function of the electron beam, this function is equal to the power absorbed by unit mass of the substance; Q pl is the energy release due to the plastic deformation per unit volume, this term is substantial for the solid parts of the target; q r and q z are the components of the heat flux in the carrying agent: where κ = κ c α c is the coefficient of heat conductivity for the mixture, κ c is the heat conductivity of the carrying agent, T c is the temperature of the carrying agent. We suppose that the volume of the dispersed phase (bubbles or droplets) remains constant in the course of the motion of the mixture as a whole and changes only due to the interfacial interaction. In this case, the volume fractions of the carrying and dispersed phases can be found from the following equations: where w is the growth rate of the dispersed phase volume per unit volume of mixture due to the interfacial interaction. Equations for the true densities of phases are as follows: where J is the growth rate of the dispersed phase mass per unit volume of mixture. Equations for internal energies of both phases take the following form: where Q is the rate of heat exchange between phases per unit volume of mixture; U tr is the internal energy of the phase, the mass of which increases in the course of phase transition: Equations (6)- (11) can be strictly derived from the conservation of volume, mass and energy in a mixture element. The functions w, J and Q express the rates of exchange between the phases. Thus, the construction of the cavitation model is reduced to the definition of these functions. One should also define pressure and stress deviators to obtain a closed system of equations. Pressure in each phase is connected with true density and internal energy through a metastable wide-range equation of state, for instance [22][23][24][25]. The model [20,21] is used to define stress deviators in the solid phase; this model takes into account the dislocation dynamics and kinetics.

Exchange rates between phases
We suppose that the dispersed phase consists of spherical particles (bubbles or droplets), which are divided into generations according to the time of nucleation. Each generation is characterized by radius R (j) , mass m (j) of one dispersed particle and concentration n (j) of particles, where the subscript (j) numerates the generations. The exchange rate of volume can be expressed as follows: where Π (j) is the change rate of the number of corresponding dispersed particles per unit volume of mixture due to nucleation or coalescence. The concentration of particles is governed by the following equation: Liquid is unstable against the liquid-vapor transition at pressure lower than the saturation vapor pressure P * v . The vapor pressure cannot be less than zero; therefore, a liquid at any negative pressure is unstable with respect to the nucleation of cavities. Nucleation of a cavity demands carrying out work against the surface tension [17]: where σ is the coefficient of surface tension; the subscript "v" relates to the vapor phase, the subscript "l" relates to the liquid phase; the record "c = l" means that the liquid is the carrying agent. The rate of nucleation can be expressed through the work W cr : where k B is the Boltzmann constant, a is an average interatomic distance, c l is the sound speed in the liquid phase; the delta symbol δ jj 0 means that the nucleation leads to arising of bubbles of the last generation j 0 with radii close to the critical radius R (j 0 ) ≈ R cr : Radii of bubbles comply with the Rayleigh-Plesset equation: where η is the melt viscosity; a dot above the symbol denotes the time derivative. Eq. (18) is valid for non-interacting bubbles and used for restricted values of volume fractions of the vapor phase: α d=v < 0.03. For higher volume fractions, another approximate equation is used: where n = j n (j) is the total concentration of cavities, α * v is the volume fraction of vapor corresponding to equilibrium with respect to pressure: where ρ * l is the density of liquid metal, at which the pressure is equal to the pressure of saturated vapor; the value ρ * l can be calculated from the equation of state. For the local volume fraction of vapor equal to one half, α v = 0.5, the percolation transition (from stage S2 to stage S3) takes place in the corresponding elements of the mixture. The following dynamics of bubbles is described by an equation similar to Eq. (18): Eq. (21) is applied at an arbitrary volume fraction of liquid droplets. The number of droplets can change as a result of coalescence in the course of the Brownian motion [26]: +0.5π where u (j) is the root-mean-square velocity of the Brownian motion: Estimation shows that the coalescence is substantial only for droplets smaller than 100 nm. The coalescence is negligible for droplets formed at ablation under the action of high-current electron irradiation.
The mass exchange rate can be expressed as follows: For stage S2 (liquid with bubbles), the kinetics of phase transition results in the following equation for the mass of a single dispersed particle (vapor bubble): where g (j) = m (j) /m 1 is the number of atoms in one bubble of generation number j, m 1 is the mass of one atom of metal, R 1 = 3 3m 1 / (4πρ * υ ), ρ * v and P * v are the density and pressure of where L j = min {0.5R j , 4λ/3}; λ is the mean free path in vapor: λ = m 1 / (Σρ v ), Σ is the effective cross-section. Eq. (26) takes into account the mode of the Knudsen flow (effusion) if the droplet size is less than the mean free path, or the continuum mode (diffusion) in the opposite case. The heat exchange between phases is described by the following equations: where L * j = min R (j) /4; 4λ/3 , C V is the specific isochoric heat capacity of vapor. Eq. (27) is for stage S2 (melt with cavities), while Eq. (28) is for stage S3 (vapor with droplets).

Parameters of the model
Surface tension is the main parameter determining the tensile strength of the melt. For little voids with a size comparable to the interatomic distance, the dependence of surface tension on radius should be taken into account, which can be approximated by Tolman's formula [14] where σ ′ is the surface tension for a flat surface. It has been shown in [27] that the surface tension of a melt linearly decreases with temperature: where σ m is the surface tension at the melting temperature T m ; the coefficient K σ is the absolute value of the derivative of surface tension over temperature, which is assumed to be a constant.
As the surface tension should be zero at the critical temperature T K , this constant can be found as K σ = σ m / (T K − T m ). According to [28], the following approximation can be used for the dependence of viscosity of liquid metals on density: where η 0 is a parameter with dimensionality of viscosity, ρ ∞ is the substance density, which corresponds to the infinitely large viscosity and zero fluidity [28]. Table 1 [27] 1358 933 T K (K) [22] 8390 8000 η 0 (µPa·s) [28] 435 91 ρ ∞ (g/cm 3 ) [28] 8.95 2.70 Figure 1. Temperature fields at sequential moments of time for the copper target. The highcurrent electron irradiation is from the left (SINUS7, R b = 1 mm); the axis of symmetry is from below. There are formation and subsequent expansion of a heated area of substance in the energy absorption zone of the electron beam.

Numerical implementation
The described model is numerically solved by the 2D version of the CRS computer program [29]. This program is designed to simulate various intensive actions on metals including high-speed impact and high-current electron irradiation. A method of division into physical processes is used with the following subproblems: i) transport of fast electrons and determination of the energy release function; ii) dynamics of the substance; iii) heat conductivity with boundary tracking of the melt and dislocation plasticity in solid parts of the metal; iv) cavitation and phase transitions in the melt. Equations for all these processes are solved independently at each time step, and the data exchange takes place at the end of each step. The transport problem for fast electrons is solved by the method [30]. The dynamics subproblem is solved by a modification of the numerical method [3]. The modification consists in eliminating the artificial viscosity and accounting for the physical viscosity instead; it allows one to obtain a stable numerical solution by using a fine enough computational grid [31]. All other equations are solved by the explicit Euler method with a varied time step. The global (hydrodynamic) time step is restricted by the Courant-Friedrichs-Lewy condition of stability. For the subproblem of phase transitions, the time step is local-it is determined separately for each computational mesh from the rates of exchange between phases; this time step can be considerably (by orders of magnitude) smaller than the global hydrodynamic step. Figure 2. Pressure fields at sequential moments of time for the copper target. The highcurrent electron irradiation is from the left (SINUS7, R b = 1 mm); the axis of symmetry is from below. One can see the formation of a high-pressure zone in the heated substance in the energy absorption zone, formation of compression and release waves and propagation of these waves towards the back surface of the target. In the layer of the substance heated by the electron irradiation, the tensile stresses (negative pressure) relax due to cavitation.

Fracture of metals in the energy absorption zone of the high-current electron beam
We numerically simulate the evolution of copper and aluminum targets under the action of a high-current electron beam with parameters corresponding to the SINUS7 facility [6]: the pulse duration at half-maximum is 50 ns, the maximum energy of electrons is 1.3 MeV, the maximum current density is 25 kA/cm 2 , and the maximum incident energy density in the center of the beam is 1.2 kJ/cm 2 . The calculations use actual experimental waveforms of accelerating voltage and current density on the target in the center of the beam. The radial distribution of current density is taken in the form of a normal distribution with an average radius R b . Calculations are performed in the 2D cylindrical formulation with normal incidence of the electron beam. Figure 3. Average density fields at sequential moments of time for the copper target. The high-current electron irradiation is from the left (SINUS7, R b = 1 mm); the axis of symmetry is from below. There are nucleation, growth and coalescence of vapor cavities in the energy absorption zone. As a result, the average density of the substance decreases here by an order of magnitude.

Evolution of the substance in the energy absorption zone
The common picture of evolution of the irradiated substance is presented in figures 1-3 for the case of copper and for the average beam radius R b = 1 mm. The target thickness is 1 mm and  Thickness of the energy absorption zone is about 0.4 mm. The substance in the energy absorption zone is heated up to 8000 K (figure 1). The heating mode is close to isochoric heating; this means that thermal expansion is negligible during the irradiation pulse. As a result, the pressure reaches the value of about 30 GPa (figure 2). Expansion of the area of high pressure leads to the formation of two compression waves, the first wave propagates towards the back surface, while the second wave runs towards the irradiated surface. Reflection of the second compression wave from the irradiated surface leads to the formation of a tensile wave with a negative value of pressure (figure 2). The action of tensile stresses on the melt heated to the temperatures of several thousands of Kelvins results in the beginning of cavitation inside the melt-nucleation of vapor bubbles; this process leads to a decrease in the average mass density of the substance (figure 3). The cavitation begins 90 ns after the onset of irradiation; the volume fraction of vapor exceeds the volume fraction of the melt in some parts of the substance beginning 180 ns after the onset of irradiation. The percolation transitions from stage S2 to stage S3 are completed in the bulk of the energy absorption zone by the time of about 300 ns. The substance beyond the energy absorption zone and a thin layer adjacent to the irradiated surface remain in the one-phase condensed state ( figure 4). As the substance in this thin surface layer is in contact with the free boundary, the tensile stresses here do not reach the value necessary for beginning of cavitation and subsequent fracture in the liquid phase.
Diameters of liquid droplets in the major part of the ablated substance inside the energy absorption zone lay in the range from several tenths of micrometers to 5 micrometers (figure 4b).
The substance at the edges of the energy absorption zone is colder, and the droplet diameters reach several tens of micrometers here.  the calculations (spatial distributions of the volume fractions of the vapor phase and the average diameter of melt droplets) for radii R b = 1 mm, 0.5 mm and 0.2 mm are shown in figures 4, 5 and 6, respectively. As follows from the presented results, the decrease in the beam radius leads to a proportional decrease in the radii of ablated areas and a change of the shape of these areas: for the beam radius of R b = 1 mm the shape of the ablated area is close to a disk with sharp edges; for the beam radius of R b = 0.2 mm the shape is closer to a sphere. The character of the spatial distribution of the droplet sizes remains the same as in the previous case: in the central part of the target, the average diameters are in the range from tenths to units of micrometers, while at the edge of the energy absorption zone the droplet diameters increase to several tens of micrometers.

Fracture of aluminum
Calculation results for aluminum are presented in figure 7. The beam parameters are the same as in the previous subsections: the pulse duration at half-maximum is 50 ns, the maximum energy of electrons is 1.3 MeV, the maximum current density is 25 kA/cm 2 , and the maximum incident energy density in the center of the beam is 1.2 kJ/cm 2 . The average radius is R b = 1 mm. The range of fast electrons in aluminum is considerably larger than in copper for the same value of the electron energy-the range is in inverse proportion to the mass density of the substance. Therefore, the electron beam heats all the thickness of the target (see figure 7a). The substance temperature reaches 3800 K in the course of irradiation. Expansion of the heated layer leads to the appearance of tensile stresses (negative pressures) with a value up to 3 GPa (see figure 7b). These tensile stresses cause active cavitation and subsequent breakup of the melt to droplets (see figures 7c and 7d). Similar to the case of copper, the formed droplets are smaller in the center of the energy absorption zone and larger near the boundary of the energy absorption zone. Most droplets have diameters in the range from 10 to 20 micrometers (figure 7d). Relatively thick layers of condensed metal with the thickness of about 0.2-0.3 mm persist along the irradiated and back surfaces.

Conclusions
A two-dimensional mathematical model of metal melt fracture at high-current electron irradiation is presented, as well as results of corresponding numerical simulations. Action of a high-current electron beam with the electron energy of about 1 MeV leads to the formation and subsequent breakup of the melt to droplets with diameters from 0.5 µm to several tens of micrometers. The finest droplets (0.5-5 µm) are generated in the center of the energy absorption zone, while the droplet sizes increase up to several tens of micrometer at the edge of the energy absorption zone, where the substance is colder. Keeping the melt in a metastable (expanded and overheated) state provides a tension wave with an amplitude of about 3-4 GPa following the shock wave and propagating deep into the target.