On experimental and numerical study of the dynamics of a liquid metal jet hit by a laser pulse

Abstract

In this paper, we present an experimental and numerical study of the laser pulse impact on the liquid metal jet target. The jet motion was recorded with the stroboscopic ultrafast shadow photography. Simulations were carried out in a two-step approach. At the first step, we simulated laser interaction with the target using the radiation hydrodynamics code 3DLINE, which accounts for a number of effects: liquid–gas phase transition, dynamics of ionization, radiation transfer, laser reflection, refraction and absorption. However, this code cannot be used on a deeply refined mesh near the liquid surface and does not account for the surface tension, which strongly affects liquid motion on the microsecond timescale after the laser pulse ends. Therefore, for the second step we employed the OpenFOAM solver, based on the volume-of-fluid method, which overcomes these limitations. The simulated target dynamics is found to be in a fairly good agreement with the experiment.

Graphic abstract

Appendices

Appendix A: 3DLINE CODE

In the current simulations, we used a one-fluid one-temperature ideal plasma model with quasi-stationary ionization. The general set of differential equations can be written in the form

$$\begin{aligned} \begin{aligned} \frac{\mathrm{d} \rho }{\mathrm{d}t}&=-\rho \nabla \cdot \vec v,\\ \rho \frac{\mathrm{d}\vec v}{\mathrm{d}t}&= -\nabla P,\\ \frac{\mathrm{d}\varepsilon }{\mathrm{d}t}&=-P \frac{\mathrm{d}(1/\rho )}{\mathrm{d}t} - \frac{1}{\rho }\nabla \cdot \vec W + G^\mathrm{rad} + G^\mathrm{las},\\ \end{aligned} \end{aligned}$$

(1)

where $\mathrm{d}/\mathrm{d}t=\partial /\partial t + \vec v\cdot \nabla$ is the material derivative, $\vec v$ is the velocity, $\rho$ is the density, $P(\rho ,T)$ is the pressure, $\varepsilon (\rho ,T)$ is the internal energy per unit mass, $\vec W=-\chi \nabla T$ is the thermoconductive heat flux, $\chi (\rho ,T)$ is the coefficient of thermal conductivity, $G^\mathrm{rad}$ is the source/sink of energy per unit mass due to the absorption/generation of thermal radiation and $G^\mathrm{las}$ is the laser energy deposition per unit mass. The thermodynamic properties were calculated as a function of density and temperature by using the FEOS code (Faik et al. 2018; Kemp and Meyer-ter Vehn 1998), based on the Thomas–Fermi model with quasi-empirical corrections in the low-temperature region.

For the radiation transport, we used the multi-group diffusion approach in quasi-stationary approximation (without the $\displaystyle \frac{1}{c}\frac{\mathrm{d}}{\mathrm{d}t}$ term which is negligible for $T\ll 1$ keV):

$$\begin{aligned} \begin{aligned} \rho G^\mathrm{rad}&=-\sum \limits _\nu \left( 4\pi j_\nu -c\varkappa ^P_{\nu } U_{\nu }\right) ,\\ \nabla \cdot \vec W_{\nu }&=4\pi j_\nu -c\varkappa ^P_{\nu } U_{\nu },\\ \vec W_{\nu }&=-\frac{c}{3\varkappa ^R_{\nu }}\nabla U_\nu . \end{aligned} \end{aligned}$$

(2)

Here, index $\nu$ denotes the spectral group, $j_\nu$ is the group emissivity, c is the speed of light, $\varkappa ^{P,R}_{\nu }$ is the mean Plank or Rosseland opacity, $U_{\nu }$ is the local radiation energy density in group $\nu$, and $\vec W_{\nu }$ is the corresponding energy flux. In this work, we used the interpolation (Novikov and Solomyannaya 1998) of the opacity and emissivity coefficients between the two tables for the limiting cases of a transparent plasma and of an optically thick plasma, both generated by the THERMOS code (Nikiforov et al. 2005).

The transport and absorption of the laser light were calculated with the hybrid model, described in detail in Refs. Basko and Tsygvintsev (2017), Tsygvintsev et al. (2016). The detailed description of the discretization scheme and the numerical algorithm can be found in Ref. Krukovskiy et al. (2017).

Appendix B: OpenFOAM PACKAGE

Unlike the 3DLINE code, the interDyMFOAM solver is based on the incompressible fluid model and does not include temperature. That is, it neglects any processes governing the temperature dynamics (such as thermal conduction, radiative heat transfer and heating by compression) as well as the temperature dependence of material properties such as its normal density, surface tension and phase transitions. At the same time, the dynamic effects of the surface tension and viscosity are fully accounted for. Gravity could also be included, but, being of negligible influence in the current setting, it was ignored in order to preserve the reflection symmetry with respect to the horizontal plane $x=0$.

Because the liquid–gas phase transition is not described by the model, the density distribution, received from the 3DLINE code, was interpreted as a mixture of two immiscible fluids, namely liquid tin and a tin vapor at low density. During the OpenFOAM simulation, these fluids are considered as one effective fluid, whose physical properties are calculated as weighted averages according to the volume-of-fluid (VOF) method (Hirt and Nichols 1981). The fluid density $\rho$ and its dynamic viscosity $\mu$ are obtained as

$$\begin{aligned} \rho = \rho _l\gamma + \rho _g(1 - \gamma ). \end{aligned}$$

$$\begin{aligned} \mu = \mu _l\gamma + \mu _g(1 - \gamma ), \end{aligned}$$

where the subscripts l and g denote the liquid and the gaseous phases, respectively. The phase fraction $\gamma$ can take values within the range $0\le \gamma \le 1$, with the values of zero and one corresponding to the regions accommodating only one phase, i.e., $\gamma =0$ for gas and $\gamma =1$ for liquid. The same approach applies to the fluid velocity

$$\begin{aligned} \vec {v} = \vec {v}_l\gamma + \vec {v}_g(1 - \gamma ). \end{aligned}$$

The evolution of the phase fraction is described by the equation

$$\begin{aligned} \frac{\partial \gamma }{\partial t} + \nabla \cdot (\vec {v}\gamma ) + \nabla \cdot [\vec {v}_r \gamma (1 - \gamma )] = 0, \end{aligned}$$

(3)

where

$$\begin{aligned} \vec {v}_r = \vec {v}_l - \vec {v}_g \end{aligned}$$

(4)

is the relative velocity vector, which should be orthogonal to the phase interface,

$$\begin{aligned} \vec {v}_r \times \nabla \gamma = 0. \end{aligned}$$

(5)

The surface tension at the liquid–gas interface generates an additional pressure gradient resulting in a force, which is evaluated per unit volume using the continuum surface force (CSF) model (Brackbill et al. 1992). For a constant surface tension $\sigma$, it can be written as

$$\begin{aligned} \vec {f}_\sigma = \sigma \kappa \nabla \gamma , \end{aligned}$$

where $\kappa$ is the mean curvature of the free surface determined from the expression

$$\begin{aligned} \kappa = - \nabla \cdot \left( \frac{\nabla \gamma }{\left| \nabla \gamma \right| }\right) . \end{aligned}$$

For incompressible fluid, the continuity equation has the form

$$\begin{aligned} \nabla \cdot \vec {v} = 0. \end{aligned}$$

(6)

Finally, with the account of everything mentioned above, the momentum equation in the Cartesian index notation takes the form

$$\begin{aligned} \begin{aligned} \frac{\partial }{\partial t}(\rho v_i) +&\frac{\partial }{\partial r_k}(\rho v_i v_k) - \frac{\partial }{\partial r_k}\left( \mu \frac{\partial v_i}{\partial r_k}\right) - \frac{\partial v_k}{\partial r_i}\frac{\partial \mu }{\partial r_k} =\\ =&-\frac{\partial p}{\partial r_i} + \sigma \kappa \frac{\partial \gamma }{\partial r_i}. \end{aligned} \end{aligned}$$

(7)

Here, p is the pressure, which is derived from the discretized form of the continuity and momentum equations (Issa 1986).

Thus, the physical model consists of the phase fraction evolution equation (3), the continuity equation (6), the momentum equation (7) and condition (5) for the “compression velocity” (4). The numerical algorithm for solving this system is based on the PISO cycle (Issa 1986), which provides the pressure–velocity coupling (Damián 2009). This coupling is achieved in the following way. At each iteration cycle, first the velocity is predicted by using the pressure field from the previous iteration. Then, by using this velocity field, the new pressure is found by solving the momentum equation, and the mass flow is recovered. At last, the velocity field is corrected according to the new mass flow. The numerical solution is performed using the finite volume method (Jasak 1996; Versteeg and Malalasekera 1995; Ferziger and Peric 2002).

Appendix C: Boundary conditions

As was mentioned before, due to the symmetry of the problem the simulations were carried out in 1/4 of the full volume, using the reflective boundary conditions in the planes $x = 0$ and $y = 0$, and the free outflow conditions at the other boundaries. For equations (1), governing matter motion and the conduction heat flow in the 3DLINE code, these conditions have the form

$$\begin{aligned} \left. P\right| _\mathrm{free\ outflow}&=0,\\ \left. \vec n\cdot \vec W\right| _\mathrm{free\ outflow}&=0,\\ \left. \vec n\cdot \vec v\right| _\mathrm{reflective}&=0,\\ \left. \vec n\cdot \vec W\right| _\mathrm{reflective}&=0,\\ \end{aligned}$$

where $\vec n$ is the unit outward normal to the bounding surface. Conductive heat flow through the domain boundaries can be neglected, because the energy is mainly carried out by the convective term $\vec v\cdot \nabla \varepsilon$ and the radiative transfer term $\vec W_\nu$. The boundary conditions for the radiation transfer equation (2) are written in the Marshack form

$$\begin{aligned} \left. \vec n\cdot \vec W_{\nu }\right| _\mathrm{reflective}&=0,\\ \left. \vec n\cdot \vec W_{\nu }\right| _\mathrm{free\ outflow}&=\frac{1}{2}c U_{\nu }.\\ \end{aligned}$$

For the laser transfer, we considered a parallel incoming beam in the $+z$ direction with the Gaussian intensity distribution

$$\begin{aligned} \left. I\right| _{z=z_{min}}=I_0(t)\cdot \exp \left( -8\frac{x^2+y^2}{D_{bw}}\right) , \end{aligned}$$

where $D_{bw}=100~\upmu$m is the diameter at the $1/e^2$ level and $I_0(t)$ is the intensity on the axis. The traced geometric rays were reflected from the $x=0$ and $y=0$ planes, and freely exited at all the other boundaries.

In the OpenFOAM, we used the most appropriate option from the available set of implemented boundary conditions. For the flow velocity $\vec {v}$, at both reflective boundaries the “zeroGradient” condition of

$$\begin{aligned} \left. \vec {v}\cdot \vec {n}\right| _\mathrm{reflective} = 0 \end{aligned}$$

(8)

was applied. For any other flow variable $\phi$ (like the phase fraction or the pressure), the “symmetryPlane” condition (Mucha 2003)

$$\begin{aligned} \left. \nabla \phi \cdot \vec {n}\right| _{y=0}= 0 \end{aligned}$$

(9)

was used in the $y=0$ plane. And in the plane $x=0$, we used the “slip” condition (Schäfer 2006)

$$\begin{aligned} \left. \vec {\tau }\cdot T \cdot \vec {n}\right| _{x=0} = 0, \end{aligned}$$

(10)

where T is the deviatoric viscous stress tensor and $\vec {\tau }$ is a unit tangential vector. In our setting, the boundary conditions (8), (9) and (8), (10) have the same physical meaning, but different numeric realization.