THEORETICAL MODEL

The HEIGHTS-PIC two-dimensional simulation package has been developed using a novel approach for the Particle-in-Cell technique with the aim of detailed analysis of intense power deposition on target materials and the resulting erosion damage, for example, over plasma facing components in Tokamaks and plasma evolution over the divertor plate. HEIGHTS-PIC is a “full-particle” representation for the momentum and energy Lagrangian variables and follows the projection of the evolution of the plasma after impingement (Genco & Hassanein, Reference Genco and Hassanein2011; Reference Genco and Hassanein2014). The position of the expanding vapor front is evaluated following each time step and compared to the original position of the birth of particles as schematically illustrated in Figure 1.

Fig. 1. (Color online) HEIGHTS-PIC modeling geometry.

One of the boundaries (the target surface) has a dynamical evolution and allows for the introduction of new generated sample particles as the material vaporizes. The sample particles or “super particles” are considered one-dimensional allowing them to be stacked along the X-axis and in layers along the Y-axis but having no dimension in this direction (Fig. 1). No map is used to transfer the properties from particles to the original grid as the conditions established by the ensemble of particles in a cell are evaluated producing the same effect inside the cell itself. The evolving characteristics of the particles inside the cells are assumed homogeneous throughout the cells. Diffusion due to the presence of a strong magnetic field is also considered and incorporated into the particle motion with possibility to include different starting diffusion coefficients and data obtained from experiments. The equations of conservation of mass, momentum, and energy are solved in two-dimensional representation (x-y plane) allowing for the plasma characteristics to be evaluated per each cell on the top of the target itself. The target heating by the laser beam is described by the equation:

(1)$${\rm \rho} c_p \displaystyle{{\partial T} \over {\partial t}} = \nabla \cdot \left({K\nabla T} \right)+ \lpar 1 - R\rpar {\rm \alpha}\, I_{S}\, e^{\lpar - {\rm \alpha} y\rpar }\comma$$

where T is the temperature, c _P is the specific heat, ρ is the density, K is the thermal conductivity, R is the target reflectivity of laser photons, α is the graphite absorption, $I_S = I_0 \exp \lpar - \vint {{\rm \kappa} dy\rpar }$ with I ₀ is maximum laser intensity, κ plasma absorption coefficient, and y is the axial coordinate.

In HEIGHTS-PIC, the solution of the transient heat transfer equations was developed with temperature and phase change dependent thermal properties (density, thermal conductivity, and specific heat). The surface temperature calculation requires the evaluation of the single portions of incident energy that go into conduction, melting, evaporation, and radiation. The kinetics of the evaporation provides the physical and mathematical connection between the out coming net particle flux from the surface and the surface temperature. Considering r as the radial distance in the plane x-z with origin in the center of the plate itself, the boundary condition at the target surface is expressed by:

(2)$$\eqalign{F\lpar r\comma \; t\rpar & = - K\lpar {T_v \lpar r \rpar } \rpar \displaystyle{{\partial T} \over {\partial z}} + {\rm \rho} \lpar {T_v \lpar r \rpar } \rpar L_v \lpar {T_v \lpar r \rpar } \rpar V\lpar {T_v \lpar r \rpar } \rpar \cr & \quad + {\rm \sigma} {\rm \varepsilon} \lpar {T_v^4 \lpar r\rpar - T_0^4 } \rpar }$$

where F(r, t) is the radial heat flux coming from the plasma, L _v is the heat of vaporization, V(T _v) is the velocity of the receding surface, and T _V = T(r, t). The velocity of the receding surface is typically a function of the instantaneous surface temperature and other materials parameters. The last term of the equation is the radiative heat transfer term. It contains T ₀, which is the temperature of the surroundings, σ is the Stefan-Boltzmann constant, and ε is material emissivity.

As the vapor starts building up and becomes thicker and thicker above the target surface, the vapor shielding of incoming laser power becomes more and more pronounced changing the boundary condition into:

(3)$$- K\displaystyle{{\partial T} \over {\partial y}}\lpar {0\comma \; t} \rpar = q_{GAS} + q_{RAD} - q_{EVAP}$$

where q _GAS is the net heat flux from the near-wall vapor cloud, q _RAD is the photon radiation heat flux absorbed at the material surface, and q _EVAP is the evaporation heat flux calculated from the enthalpy of evaporation for the material under investigation, i.e., carbon in this case. The evaporation flux is calculated using the Hertz-Knudsen-Langmuir theory of evaporation and condensation under non-equilibrium conditions (Hassanein et al., Reference Hassanein, Kulcinski and Wolfer1981; Reference Hassanein, Kulcinski and Wolfer1982). As the laser pulse comes to an end, the plasma starts cooling down and starts its expansion phase. The expansion of the vapor/plasma into vacuum is then determined by solving the MHD equations for conservation of mass, momentum, and energy (4a) associated with Maxwell equations (4b):

(4a)$$\eqalign{& \displaystyle{{\partial {\rm \rho} } \over {\partial t}} + \mathop \nabla \limits^ \to {\rm \rho} \vec V = 0 \cr & \displaystyle{{\partial {\rm \rho} \vec V_x } \over {\partial t}} + \mathop \nabla \limits^ \to {\rm \rho} V_x \vec V + \displaystyle{\partial \over {\partial x}}\left({P + \displaystyle{{B^2 } \over {2{\rm \mu} _0 }}} \right)- \displaystyle{1 \over {{\rm \mu} _0 }}\left({\mathop {B}\limits^ {\to}} {\mathop {\nabla} \limits^{\to} } \right)B_x = F_x \cr & \displaystyle{{\partial {\rm \rho} \vec V_y } \over {\partial t}} + \mathop \nabla \limits^ \to {\rm \rho} V_y \vec V + \displaystyle{\partial \over {\partial x}}\left({P + \displaystyle{{B^2 } \over {2{\rm \mu} _0 }}} \right)- \displaystyle{1 \over {{\rm \mu} _0 }}\left({\mathop {B}\limits^{\to} \mathop {\nabla} \limits^{\to} } \right)B_y = F_y = 0 \cr & \displaystyle{{\partial {\rm \rho} \vec V_z } \over {\partial t}} + \mathop \nabla \limits^ \to {\rm \rho} V_z \vec V - \displaystyle{1 \over {{\rm \mu} _0 }}\left({\mathop {B}\limits^{\to}} {\mathop {\nabla} \limits^{\to} } \right)B_z = 0 \cr & \displaystyle{\partial \over {\partial t}}\left\{{{\rm \rho} e_i + \displaystyle{{{\rm \rho} V^2 } \over 2} + \displaystyle{{B^2 } \over {2{\rm \mu} _0 }}} \right\}+ \mathop \nabla \limits^ \to \left\{\matrix{\vec V_{} \left({{\rm \rho} e_i + \displaystyle{{{\rm \rho} V^2 } \over 2}} \right)+ \hfill \cr + \,{\mathop{P}\limits^{\to}} \cdot \vec V - K\mathop \nabla \limits^ \to T_e + S_{Rad} + \displaystyle{{\mathop E\limits^ \to \times \mathop B\limits^ \to } \over {{\rm \mu} _0 }} \hfill} \right\}\cr &\quad= Q_{Beam}}$$

(4b)$$\eqalign{& \displaystyle{{\partial \mathop B\limits^ \to } \over {\partial t}} = - \mathop \nabla \limits^ \to \times \mathop E\limits^ \to \cr & \displaystyle{1 \over {{\rm \mu} _0 }}\mathop \nabla \limits^ \to \times \mathop B\limits^ \to = \mathop J\limits^ \to \cr & \mathop \nabla \limits^ \to \bullet \mathop B\limits^ \to = 0 \cr & \mathop \nabla \limits^ \to \bullet \mathop E\limits^ \to = {\rm \rho} _e \cr & \mathop E\limits^ \to + \mathop V\limits^ \to \times \mathop B\limits^ \to = \displaystyle{1 \over {\rm \sigma} }\mathop J\limits^ \to }$$

where ρ is the density, $\vec V$ is the vapor velocity, P is the pressure, E is the energy, K is the vapor conductivity, Q _b is the incident photons flux from the incoming laser pulse, and Q _r is the radiation flux. The equation describing the vapor motion in a magnetic field is given by:

(5)$${\rm \rho} \displaystyle{{\partial \mathop V\limits^ \to } \over {\partial t}} = - \nabla P + \vec J \times \vec B$$

where $\vec J \times \vec B$ is the magnetic force. The vapor plasma, once it is ionized, moves freely along the magnetic field lines. More details of the model can be found in previous publications (Hassanein, Reference Hassanein1994; Hassanein & Konkashbaev, Reference Hassanein and Konkashbaev1995).

NUMERICAL SIMULATION

The performance of a graphite target irradiation was simulated numerically mimicking a ND:YAG laser operating with a Gaussian beam shape at it fundamental wavelength of 1064 nm with pulse energy of 400 mJ, total pulse duration of 30 ns with full width at half maximum of 10 ns. The laser spot size on the target has a diameter of 0.1784 cm (area of 2.5 mm²) with a laser total energy flux close to 16 J/cm² with intensity of 1.6 GW/cm². Plasma evolution was simulated up to 40 ns from the start of the laser beam. A magnetic field of 0.5 T was assumed perpendicular to the target. The number of super particles initially loaded into HEIGHTS-PIC was relatively low (N = 20) and the numerical grid used was made of 6 cells along the X-axis and 10 cells along the Y-axis for a total cell dimension of 0.0296 cm × 0.033 cm.

RESULTS AND DISCUSSION

The effect of the imping laser onto the target surface and plasma energy absorption is described in Figures 2 and 3. The temperature of the target surface reaches its maximum around 8060 K at the center of the beam (r = 0) in the x-z plane and is lower for increasing values of the radius r, which follows the laser beam profile. This maximum temperature value is very close to recent work for similar laser energy flux (Bulgakova et al., Reference Bulgakova, Bulgakov and Babich2004).

Fig. 2. (Color online) Evolution of the calculated target temperature along time at different locations.

Fig. 3. (Color online) Power fluxes in the first instants of the laser pulse.

The effective power flux reaching the surface changes rapidly with time. Given the temporal Gaussian evolution of the intensity of the laser beam, it is evident that initially all the energy goes into the target plate increasing its temperature. The value of laser reflectivity used is taken from (Bulgakova et al., Reference Bulgakova and Bulgakov2001; Malvezzi & Bloembergen, Reference Malvezzi and Bloembergen1986) and provides the variation of carbon reflectivity as a function of the surface temperature. It is equal to R = 0.21 at the beginning (T = 300 K) and becomes very small at higher temperatures (less than 0.05 for T > 6000 K). When plasma starts being formed around 11 ns, the carbon plate has already reached a surface temperature close to 5000 K so that its reflectivity according to Bulgakova and Bulgakov (Reference Bulgakova and Bulgakov2001) is calculated roughly equal to 0.08 and rapidly diminishing.

In fact, in the very first nanoseconds of the pulse, and up to the time in which the temperature in the target has been raised enough to vaporize the material and eject particles from the surface (i.e., 10.8 ns), the total laser power flux goes into the plate and it is fully absorbed (Fig. 3). As carbon plasma is formed and vapor particles are produced and emitted into the system, some of the incoming laser power flux starts being absorbed into the plasma itself through photon interaction, reducing the effective amount of power reaching the plate. As this phenomenon becomes more and more important, more and more energy is absorbed into the plume effectively screening the target for further damage. The direct laser intensity to target drops rapidly toward zero up to a point in which the target is totally screened from it and further erosion is eventually produced by the radiation only as shown in Figure 4.

Fig. 4. (Color online) Calculated evolution of the erosion of the carbon target center along time and variation of maximum T surface at the center of the plate.

The value of the maximum erosion calculated after one laser pulse (0.33 µm) is very close to the ones provided in the literature both theoretically and experimentally for similar conditions (Bulgakov & Bulgakova, Reference Bulgakov and Bulgakova1999; Hoffman et al., Reference Hoffman, Moscicki and Szymanski2011). It is also very similar to the value of 0.40 µm predicted by the well-known and well benchmarked HEIGHTS-LPP package (Sizyuk et al., Reference Sizyuk, Hassanein and Sizyuk2006).

Looking at the integrated energy fluxes at the end of the pulse (Fig. 5), it can be seen that of the total energy provided (100%, blue line) only 22% is effectively absorbed into the plate while all the rest (78% black line) is absorbed into plasma (45%, cyan line) or reradiated back (33%, green line) This distribution of total energies is also very close to the one provided by other theoretical models (Bulgakova et al., Reference Bulgakova, Bulgakov and Babich2004).

Fig. 5. (Color online) Calculated integrated energy fluxes of the system.

Values of the electron density and temperature of the plasma are calculated along the entire mesh and for each time step. At 40 ns the profile of temperature close to the target still has a Gaussian shape. The peak of temperature is around 6.6 eV and drops with distance as expected (Harilal et al., Reference Harilal, Bindhu, Issac, Nampoori and Vallabhan1997). The comparison with experimental data for similar conditions is quite good given the experimental error of ±15% for the electron temperature (Fig. 6).

Fig. 6. (Color online) Calculated plasma temperature profiles close to the target at 40 ns and comparison with experimental data (Hoffman et al., Reference Hoffman, Moscicki and Szymanski2011).

Electron density has a similar profile. The experimental values provided in the literature report at such close distance to target an error of the measurements around 50–60% for electron density (Hoffman et al., Reference Hoffman, Moscicki and Szymanski2011) so that the comparison can be considered fair (Fig. 7).

Fig. 7. (Color online) Calculated electron density profiles close to the target at 40 ns and comparison with experimental data (Hoffman et al., Reference Hoffman, Moscicki and Szymanski2011).

It should be pointed out, that HEIGHTS-PIC has, into the solution of the hydrodynamics equations, the implementation of the contribution of the magnetic force. As observed already (Pathak & Chandy, Reference Pathak and Chandy2009; Neogi & Thareja, Reference Neogi and Thareja2001) the presence of the $\vec J \times \vec B$ magnetic force changes the behavior of the plasma plume with respect to the one expanding for example into vacuum having a strong influence on the ionization characteristics as well as on the velocity of the plume itself. This translates for axial magnetic field as for HEIGHTS-PIC in a stronger confinement close to target.

Since some uncertainty still exists regarding the reflectivity of graphite exposed to laser pulse, Tokamak graphite has also been considered. The reflectivity of dense graphite has been estimated at room temperature around R = 0.40 to 0.29 while it appears that due to high porosity R = 0.22 (Semerok et al., Reference Semerok, Fomichev, Weulersse, Brygo, Thro and Grisolia2007) in Tokamaks; this value is also almost independent from surface temperature (Semerok et al., Reference Semerok, Fomichev, Weulersse, Brygo, Thro and Grisolia2007). For this reason, HEIGHTS-PIC was also run with R = 0.22 constant and R = 0.40 for the entire simulation in order to evaluate differences with the previous examined case. Particularly interesting are the profiles of erosion. In fact as shown in Figure 8 final erosions can be different according to which type of graphite reflectivity is used into the calculation.

Fig. 8. (Color online) Calculated evolution of the erosion of the carbon target center along time with different values of the carbon reflectivity.

With variable reflectivity the calculated erosion is more significant than the other two cases with constant reflectivity. In fact as soon as the temperature is elevated enough from the original 300 K, R tends to go down and becomes almost insignificant at high temperatures. Instead for the other two cases with R constant, since more energy is reflected back and eventually lost, less is the energy reaching the surface producing less erosion. So a more accurate evaluation of the behavior of reflectivity along with time under beam irradiation is critical in order to establish more accurate erosion and lifetime values.