Fluid PIC scheme

We use the Harlow’s splitting PIC scheme to solve the hydrodynamic equations. Without loss of generality, the equations can be written in the abstract operator form: $\frac{\partial q}{\partial t}+\widehat{A}q=\mathrm{0,}$ (1) where q(r,t) is an arbitrary hydrodynamic quantity and $\widehat{A}$ is an abstract operator. According to the rule of operator decomposition, the solution of equation (1) at time step $\text{Δ}t$ is reduced to a sequential solution of the two auxiliary problems [28] $\{\begin{array}{l}\frac{\partial \tilde{q}\left(r\mathrm{,}t\right)}{\partial t}+\widehat{E}\tilde{q}\left(r\mathrm{,}t\right)=\mathrm{0,}\hfill \\ \frac{\partial q\left(r\mathrm{,}t\right)}{\partial t}+\widehat{L}q\left(r\mathrm{,}t\right)=\mathrm{0,}\hfill \end{array}$ (2) where $\widehat{A}=\widehat{E}+\widehat{L}$. The first equation in Eq. 2 corresponds to the “Euler step” and the second corresponds to “Lagrange step” in the Harlow’s splitting scheme, respectively. In the “Euler step”, the operator $\widehat{E}$ does not incorporate the spatial divergence operator. Thus, the equation is easily solved on a fixed spatial grid. In the “Lagrange step”, pseudo-particles are introduced to carry mass density, momentum density, and specific internal energy density. In one time step, as pseudo-particles move to new positions, new hydrodynamic quantities can be obtained on the grids by summing the pseudo-particles. The “Lagrange step” can be seen as a computational procedure for modeling particle migration, which compensates for the transport effect that is neglected in the “Euler step”.

Without loss of generality, the second equation in Eq. 2 can be written as follows: $\frac{\partial q}{\partial t}+\nabla \cdot (qU)=\mathrm{0,}$ (3) where U is the flow velocity, and q=(p,pU, pe) represent the mass, momentum and internal energy density, respectively. Eq. 3 is in the form of the conservation equation, the solution of which can be written as a summation of pseudo-particles. $q(r\mathrm{,}t)={\displaystyle \sum _{j=1}^N}{Q}_{j}R(r\mathrm{,}{r}_j(t))\mathrm{.}$ (4) The total number of pseudo-particles is denoted by N, and Q_j represents the hydrodynamic quantities carried by the jth pseudo-particle. The value of Q_j remains constant in one “Lagrange step”. The kernel function of the pseudo-particle denoted as R(r,r_j(t)), depends on the current space coordinate r and the radius vector r_j of the jth pseudo-particle center. This kernel function satisfies certain universal properties of typical $\{\begin{array}{l}R\left(r_1\mathrm{,}{r}_2\right)=R\left(r_2\mathrm{,}{r}_1\right)⩾\mathrm{0,}\hfill \\ \frac{\partial R}{\partial r_1}=-\frac{\partial R}{\partial r_2}\mathrm{,}\hfill \\ {\displaystyle {\int }_{\text{Ω}}}\text{d}{r}_{1}R\left(r_1\mathrm{,}{r}_2\right)=\mathrm{1,}\hfill \end{array}$ (5) where Ω denotes the full space. Given the aforementioned characteristics and any smooth finite function, the representation q(r,t) in Eq. 4 simplifies Eq. 3 to satisfy the equations of motion for pseudo-particles, $\frac{\text{d}{r}_j}{\text{d}t}=U(r_j(t))\mathrm{,}j=\mathrm{1,2,}\dots \mathrm{,}N\mathrm{.}$ (6) Our code comprises three primary modules: the energy deposition module, equation of state (EOS) module, and ideal hydrodynamic module, which are complemented by a post-processing script for the blow-off impulse calculation to constitute the complete Xablation2D code.

Energy deposition

The energy deposition module is responsible for estimating the energy transfer between the X-rays and matter. At a microscopic level, X-rays primarily interact with matter through electrons. Initially, the energy of the photons is transferred directly to the electrons, which then interact with the atoms in the matter, resulting in energy deposition. The photoelectric effect is dominant in the low-energy region. As the photon energy increases, Compton scattering contributes to the energy deposition [29]. In far-field X-ray irradiation, the degree of ionization of the material is so low that the effects of the plasma can be disregarded. The energy flux for parallel X-rays incident on the target material can be estimated by considering the distance x traveled in the direction of incidence. $\text{Φ}(x)={\text{Φ}}_0\frac{{\displaystyle {\int }_0^{\infty }}f\left(\lambda \mathrm{,}T\right)\exp \left[-\mu (\lambda )\rho x\right]\text{d}\lambda }{{\displaystyle {\int }_0^{\infty }}f\left(\lambda \mathrm{,}T\right)\text{d}\lambda }\mathrm{,}$ (7) where Φ₀ and Φ(x) are energy flux at the initial position and at position x, respectively; $f(\lambda \mathrm{,}T)$ is the energy spectrum of X-rays, which depends on the photon wavelength λ and radiation temperature T; ρ is the mass density; $\mu (\lambda )$ is the mass absorption coefficient associated with the photon wavelength. The energy spectrum $f(\lambda \mathrm{,}T)$ approximates a blackbody spectrum for X-rays produced by a nuclear explosion. In numerical simulations, it is necessary to truncate and discretize the energy spectrum. Hence, the expression for the energy flux in computing can be written as $\text{Φ}(x)={\text{Φ}}_0{\displaystyle \sum _j{w}_j}\exp \left[-\mu \left({\lambda }_j\right)\rho x\right]\mathrm{,}$ (8) where the subscript j is the discretized energy group index and wj represents the proportion of the incident energy flux of monochromatic light with a specific wavelength λj. The quantity eR signifies the amount of photon energy deposited through the interaction between X-rays and matter per unit mass and per unit time within the region $x\sim x+\text{Δ}x$ $e_R=\frac{\text{Φ}(x)-\text{Φ}(x+\text{Δ}x)}{\rho \text{Δ}x\tau },$ (9) where τ denotes a time increment. It is postulated in the subsequent discussion that all the deposited energy is converted into internal energy.

Equation of state

X-ray irradiation induces significant changes in the state of matter, including phase transitions, thermal expansion, and shock compression, which results in a large range of material parameters. Consequently, it is necessary to employ distinct equations of state to characterize the expansion and compression regions of a material. For the thermal expansion region, the PUFF EOS [15] was adopted. $\begin{array}{l}p=\rho \left[\gamma -1+\left({\Gamma }_0-\gamma +1\right)\sqrt{\frac{\rho }{{\rho }_0}}\right]\\ \left[e-e_{\text{s}}\left\{1-\exp \left[\frac{N{\rho }_0}{\rho }\left(1-\frac{{\rho }_0}{\rho }\right)\right]\right\}\right]\mathrm{,}\rho <{\rho }_0\mathrm{,}\end{array}$ (10) where ρ₀ and ρ are the initial density and current mass densities, respectively, p is the material pressure, Γ₀ is the Güneisen coefficient, γ is the specific heat ratio of vaporized gas, e_s is the sublimation energy, and $N=C_0^2/{\Gamma }_0{e}_{\text{s}}\cdot C_0$ is a Hugoniot parameter that determines the shock wave velocity D in a solid material with postshock velocity u and another Hugoniot parameter λ. The relation is $D=C_0+\lambda u\mathrm{.}$ (11) For the compression zone, the Güneisen EOS on the Hugoniot line is used, which is expressed as [30] $p=p_{\text{H}}\left(v\right)+{\rho }_0{\Gamma }_0\left(e-e_{\text{H}}\right)\mathrm{,}\rho ⩾{\rho }_0\mathrm{,}$ (12) where $p_{\text{H}}\left(v\right)$ and e_H are $p_{\text{H}}\left(v\right)=\frac{{\rho }_0{C}_0^2\left(1-v/v_0\right)}{\left[1-\lambda {\left(1-v/v_0\right)}^2\right]}$ (13) and $e_{\text{H}}=\frac{1}{2}{p}_{\text{H}}\left(v_0-v\right)\mathrm{,}$ (14) respectively, representing the postshock pressure and specific internal energy when the preshock is stationary; $v=1/\rho$ is the specific volume and v₀ is the initial specific volume.

In addition to the analytic expression of the EOS, an open-source code, FEOS, provides the EOS for a wide range of temperatures and densities in tabular form [31]. The FEOS, based on the quotidian equation of state (QEOS) model [32], calculates the material-specific Helmholtz free energy $F(\rho \mathrm{,}T)$ directly, which is widely used in computational fluid dynamics codes.

Ideal hydrodynamics

In the hydrodynamic module, the effects of thermal radiation and heat conduction can be ignored. The far-field X-ray flux with respect to these issues is $\mathcal{O}(100{\text{ J/cm}}^2)$. In this case, the material temperature is $\mathcal{O}(1\text{ eV})$ and the corresponding radiation pressure is only $\mathcal{O}(1\text{Pa})$, which is far lower than the material pressure. The velocity of the thermal shock wave in solid materials is approximately $\mathcal{O}(1\text{km/s})$, which is much faster than the heat conduction; therefore, these effects can be ignored. Furthermore, our code aims to calculate the blow-off impulse occurring in the thermal expansion-vaporized region. This region is characterized by a significantly lower material stress than the material pressure; thus, we can also ignore the stress for the impulse calculation in the following.

The 2D governing equations for ideal hydrodynamics are as follows: $\{\begin{array}{l}\frac{\partial \rho }{\partial t}+\frac{\partial }{\partial x}\left(\rho u_x\right)+\frac{\partial }{\partial y}\left(\rho u_y\right)=\mathrm{0,}\hfill \\ \frac{\partial }{\partial t}\left(\rho u_x\right)+\frac{\partial }{\partial x}\left(\rho u_x^2\right)+\frac{\partial }{\partial y}\left(\rho u_x{u}_y\right)=-\frac{\partial p}{\partial x}\mathrm{,}\hfill \\ \frac{\partial }{\partial t}\left(\rho u_y\right)+\frac{\partial }{\partial y}\left(\rho u_y^2\right)+\frac{\partial }{\partial x}\left(\rho u_x{u}_y\right)=-\frac{\partial p}{\partial y}\mathrm{,}\hfill \\ \frac{\partial w}{\partial t}+\frac{\partial }{\partial x}\left(w{u}_x\right)+\frac{\partial }{\partial y}\left(w{u}_y\right)=e_R\mathrm{,}\hfill \end{array}$ (15) where ρ is the mass density, $u=(u_x\mathrm{,}{u}_y)$ is the fluid velocity, p is the pressure, e is the specific internal energy, $w=\rho \left(e+u^2/2\right)$ is the total energy, and eR is the energy deposited by X-ray irradiation per unit mass per unit time. According to Harlow’s splitting scheme, the 2D governing equations can be divided into the following distinct groups: $\{\begin{array}{l}\frac{\partial {\rho }_1}{\partial t}=\mathrm{0,}\hfill \\ \frac{\partial }{\partial t}\left({\rho }_1{u}_{x\mathrm{,1}}\right)=-\frac{\partial p}{\partial x}\mathrm{,}\hfill \\ \frac{\partial }{\partial t}\left({\rho }_1{u}_{y\mathrm{,1}}\right)=-\frac{\partial p}{\partial y}\mathrm{,}\hfill \\ \frac{\partial w}{\partial t}=-\frac{\partial }{\partial x}\left(p{u}_{x\mathrm{,1}}\right)-\frac{\partial }{\partial y}\left(p{u}_{y\mathrm{,1}}\right)+e_R\mathrm{,}\hfill \end{array}$ (16) and $\{\begin{array}{l}\frac{\partial {\rho }_2}{\partial t}+\frac{\partial }{\partial x}\left({\rho }_2{u}_{x\mathrm{,2}}\right)+\frac{\partial }{\partial y}\left({\rho }_2{u}_{y\mathrm{,2}}\right)=\mathrm{0,}\hfill \\ \frac{\partial }{\partial t}\left({\rho }_2{u}_{x\mathrm{,2}}\right)+\frac{\partial }{\partial x}\left({\rho }_2{u}_{x\mathrm{,2}}^2\right)+\frac{\partial }{\partial y}\left({\rho }_2{u}_{x\mathrm{,2}}{u}_{y\mathrm{,2}}\right)=\mathrm{0,}\hfill \\ \frac{\partial }{\partial t}\left({\rho }_2{u}_{y\mathrm{,2}}\right)+\frac{\partial }{\partial x}\left({\rho }_2{u}_{x\mathrm{,2}}{u}_{y\mathrm{,2}}\right)+\frac{\partial }{\partial y}\left({\rho }_2{u}_{y\mathrm{,2}}^2\right)=\mathrm{0,}\hfill \\ \frac{\partial w_2}{\partial t}+\frac{\partial }{\partial x}\left(u_{x\mathrm{,2}}{w}_2\right)+\frac{\partial }{\partial y}\left(u_{y\mathrm{,2}}{w}_2\right)=\mathrm{0,}\hfill \end{array}$ (17) where subscripts 1 and 2 denote “Euler step” and “Lagrange step”, respectively.

Blow-off impulse calculation

When X-rays irradiate the surface of a material, the material undergoes sublimation, forming an evaporation layer. The vaporized material is violently ejected outward into the surroundings, generating a blow-off impulse. According to the impulse theorem, the specific impulse in a given direction is equal to the momentum change $I={\displaystyle {\int }_{P_0}^P}\text{d}P=P-P_0=m\left(u-u_0\right)\mathrm{,}$ (18) where $P\mathrm{,}u\mathrm{,}{P}_0$ and u₀ represent the final momentum, final flow velocity, initial momentum, and initial flow velocity, respectively, in a certain direction, and m is the material mass. Typically, the initial velocity of the vaporized material is zero, that is u₀=0. For the numerical calculations, Eq. 18 can be written in summation form [33] $I={\displaystyle \sum _{u_j\langle \mathrm{0,}{e}_j\rangle e_{\text{s}}}{m}_j\left|u_j\right|}\mathrm{,}$ (19) where mj and $u_j$ are the material mass and velocity in the jth grid, respectively, and accounts is the sum of the momentum in all grids in which the vaporized material ejects outward. However, BOI can also be computed using the definition: $I={\displaystyle {\int }_{t_0}^t}{p}_{\text{g}}\text{d}t\mathrm{,}$ (20) where p_g represents the pressure exerted by the ejected gas on the surface of the solid material at rest and t and t₀ are the final and initial times, respectively. For the numerical calculations, the discrete form of Eq. 20 is as follows: $I={\displaystyle \sum _n{p}_{\text{g}}\text{Δ}t}\mathrm{.}$ (21) where n represents the nth vaporized grid.

Initialization

We discretize the simulation domain into rectangular grids in Cartesian coordinates. The physical quantities of the fluid, except for the mass density, are initialized by manually assigning them to grid points. The mass density is initialized by summing the weights of the pseudo-particles, as follows: $\begin{array}{c}{\rho }_{i\mathrm{,}j}={\rho }_{i\mathrm{,}j}+{\rho }_p\cdot {\omega }_{i\mathrm{,}j}\mathrm{,}\\ {\rho }_{i+\mathrm{1,}j}={\rho }_{i+\mathrm{1,}j}+{\rho }_p\cdot {\omega }_{i\mathrm{,}j}\mathrm{,}\\ {\rho }_{i\mathrm{,}j+1}={\rho }_{i\mathrm{,}j+1}+{\rho }_p\cdot {\omega }_{i\mathrm{,}j}\mathrm{,}\\ {\rho }_{i+\mathrm{1,}j+1}={\rho }_{i+\mathrm{1,}j+1}+{\rho }_p\cdot {\omega }_{i\mathrm{,}j}\mathrm{,}\end{array}$ (22) where ${\rho }_{i\mathrm{,}j}$ denotes the mass density of the grid, ρp denotes the mass density of the pseudo-particle, and ωi,j denotes the weight factor. Notably, in our code, the value of ρp remains constant throughout the simulation. The area-weighting method illustrated in Fig. 1 is employed to distribute the mass density in the 2D simulation, and weight factors are expressed as $\begin{array}{c}{\omega }_{i\mathrm{,}j}=\frac{A_1}{A_1+A_2+A_3+A_4}\mathrm{,}\\ {\omega }_{i+\mathrm{1,}j}=\frac{A_2}{A_1+A_2+A_3+A_4}\mathrm{,}\\ {\omega }_{i\mathrm{,}j+1}=\frac{A_3}{A_1+A_2+A_3+A_4}\mathrm{,}\\ {\omega }_{i+\mathrm{1,}j+1}=\frac{A_4}{A_1+A_2+A_3+A_4}\mathrm{,}\end{array}$ (23) where $A_i\left(i=\mathrm{1,2,3,4}\right)$ represents the area of overlap between the pseudo-particle and the grid.

Fig. 1Full-screen View

Area-weighting method for the pseudo-particle in 2D simulation.The circle p represents the pseudo-particle center, and the black dot (i,j) is the grid point

Euler step

The radiative deposition energy and the equation of state are discretized on the grid as follows: $e_{R\mathrm{,}i\mathrm{,}j}^n=\frac{{\text{Φ}}_{i-\mathrm{1,}j-1}^n-{\text{Φ}}_{i\mathrm{,}j}^n}{{\rho }_{i\mathrm{,}j}^n{h}_1\tau }\mathrm{,}$ (24) and $(\begin{array}{l}\begin{array}{c}{p}_{i\mathrm{,}j}^n=p_{H\mathrm{,}i\mathrm{,}j}^n\\ \begin{array}{cc}+{\Gamma }_0{\rho }_0\left[e_{i\mathrm{,}j}^n-\frac{1}{2}{p}_{H\mathrm{,}i\mathrm{,}j}^n\left(v_0-v_{i\mathrm{,}j}^n\right)\right]\mathrm{,}& {\rho }_{i\mathrm{,}j}⩾{\rho }_0\mathrm{,}\end{array}\end{array}\hfill \\ \begin{array}{c}{p}_{i\mathrm{,}j}^n={\rho }_{i\mathrm{,}j}^n\\ \begin{array}{cc}\times \left[\gamma -1\left({\Gamma }_0-\gamma +1\right)\sqrt{\frac{{\rho }_{i\mathrm{,}j}^n}{{\rho }_0}}\right]e_{A\mathrm{,}i\mathrm{,}j}^n\mathrm{,}& {\rho }_{i\mathrm{,}j}<{\rho }_0\mathrm{,}\end{array}\end{array}\hfill \end{array}$ (25) where $p_{H\mathrm{,}i\mathrm{,}j}^n$ and $e_{A\mathrm{,}i\mathrm{,}j}^n$ are: $p_{H\mathrm{,}i\mathrm{,}j}^n=\frac{{\rho }_0{C}_0^2\left(1-v_{i\mathrm{,}j}^n/v_0\right)}{\left[1-\lambda {\left(1-v_{i\mathrm{,}j}^n/v_0\right)}^2\right]}\mathrm{,}$ (26) and $e_{A\mathrm{,}i\mathrm{,}j}^n=e_{i\mathrm{,}j}^n-e_{\text{s}}\left\{1-\exp \left[\frac{N{\rho }_0}{{\rho }_{i\mathrm{,}j}^n}\left(1-\frac{{\rho }_0}{{\rho }_{i\mathrm{,}j}^n}\right)\right]\right\}\mathrm{,}$ (27) respectively. Subsequently, the hydrodynamic equations (Eq. 16) in the “Euler step” can be discretized as $\{\begin{array}{l}{u}_{x\mathrm{,}i\mathrm{,}j}^{n+1}=u_{x\mathrm{,}i\mathrm{,}j}^n-\frac{\tau }{h_1{\rho }_{i\mathrm{,}j}^n}\left(p_{i+1/\mathrm{2,}j}^n-p_{i-1/\mathrm{2,}j}^n\right)\mathrm{,}\hfill \\ u_{y\mathrm{,}i\mathrm{,}j}^{n+1}=u_{y\mathrm{,}i\mathrm{,}j}^n-\frac{\tau }{h_2{\rho }_{i\mathrm{,}j}^n}\left(p_{i\mathrm{,}j+1/2}^n-p_{i\mathrm{,}j-1/\mathrm{2,}j}^n\right)\mathrm{,}\hfill \\ \begin{array}{c}{e}_{i\mathrm{,}j}^{n+1}=e_{i\mathrm{,}j}^n-\frac{\tau p_{i\mathrm{,}j}^n}{{\rho }_{i\mathrm{,}j}^n}(\frac{u_{x\mathrm{,}i+1/\mathrm{2,}j}^n-u_{x\mathrm{,}i-1/\mathrm{2,}j}^n}{h_1}\\ +\frac{u_{y\mathrm{,}i\mathrm{,}j+1/2}^n-u_{y\mathrm{,}i\mathrm{,}j-1/2}^n}{h_2})+e_{R\mathrm{,}i\mathrm{,}j}^n\mathrm{,}\end{array}\hfill \end{array}$ (28) where the subscript (i,j) denotes the grid index in the spatial coordinate system, the superscript n represents the nth time step, τ is the time increment, h₁ and h₂ are the grid sizes in the x and y directions, respectively. The mass density remained constant, as shown in Eq. 16, that is, ${\rho }_{i\mathrm{,}j}^{n+1}={\rho }_{i\mathrm{,}j}^n$.

Lagrange step

In the “Lagrange step”, the fluid is divided into many pseudo-particles with finite sizes. We assume an internal distribution function of the fluid quantities within a pseudo-particle [34] $\begin{array}{c}q\left(\xi \mathrm{,}\eta \right)=\left\{\left(q_{i+\mathrm{1,}j+1}-q_{i+\mathrm{1,}j}\right)\left(2\eta -\delta \eta \right)+q_{i+\mathrm{1,}j}\right\}\left(2\xi -\delta \xi \right)\\ -\left\{\left(q_{i\mathrm{,}j+1}-q_{i\mathrm{,}j}\right)\left(2\eta -\delta \eta \right)+q_{i\mathrm{,}j}\right\}\left(2\xi -\delta \xi -1\right)\mathrm{,}\end{array}$ (29) where the pseudo-particle size is normalized and q(ξ,η) is the physical quantity carried by the pseudo-particles¹, the coordinates (ξ,η) represent the internal position within the pseudo-particle with the origin located at the bottom-left corner and (δξ, δη) are intervals from the origin of the internal coordinate system to the boundaries of the grid in the x and y directions, respectively. The results are shown in Fig. 2, where the red lines represent the grids, the red dot represents the pseudo-particle center, and the blue square represents the pseudo-particle size.

Fig. 2Full-screen View

Distribution of quantities within a pseudo-particle. The red dot p represents the pseudo-particle center, the blue square represents the pseudo-particle size, and the black dot is the grid point. The label ^* represents the interval when the pseudo-particle reaches a new position

q_p represents the physical quantity carried by the pseudo-particle located at the center of the internal coordinates, which can be obtained by integrating q(ξ,η) from 0 to 1. $\begin{array}{c}{q}_{\text{p}}={\displaystyle \underset{0}{\overset{1}{\int }}{\displaystyle \underset{0}{\overset{1}{\int }}q\left(\xi \mathrm{,}\eta \right)}}\text{d}\xi \text{d}\eta \\ =q_{i+\mathrm{1,}j+1}\left(1-\delta \xi \right)\left(1-\delta \eta \right)\\ +q_{i+\mathrm{1,}j}\left(1-\delta \xi \right)\delta \eta +q_{i\mathrm{,}j+1}\delta \xi \left(1-\delta \eta \right)+q_{i\mathrm{,}j}\delta \xi \delta \eta \\ =q\left(\frac{1}{2}\mathrm{,}\frac{1}{2}\right)\mathrm{.}\end{array}$ (30) The equations for pseudo-particle motion are $\{\begin{array}{ll}{x}_p^{n+1}\hfill & =x_p^n+\tau u_{p\mathrm{,}x}\mathrm{,}\hfill \\ y_p^{n+1}\hfill & =y_p^n+\tau u_{p\mathrm{,}y}\mathrm{,}\hfill \end{array}$ (31) where ($x_p\mathrm{,}{y}_p$) denotes the position of the pseudo-particle center and $\left(u_{p\mathrm{,}x}\mathrm{,}{u}_{p\mathrm{,}y}\right)$ represents the velocity of the pseudo-particle center, as determined by Eq. 30. When the pseudo-particle reaches a new position as shown in Fig. 2, we recompute the new intervals from the origin of the internal coordinate, denoted as $\left(\delta {\xi }^{*}\mathrm{,}\delta {\eta }^{*}\right)$. By summing the pseudo-particles at new positions, new fluid quantities can be derived on the grid. In this study, we present three algorithms for the summation process: the area-weighting method (AWM), integration-weighting method (IWM), and interpolation–integration–weighting method (IIWM), each of which exhibits different types of numerical diffusion [34]. The details of summation algorithms are shown in Appendix 6. The three algorithms have different numerical diffusions and noise and can be selected according to the requirements of the actual problems.

Parallelization

The message-passing interface (MPI) is employed in the construction of the parallelization [35]. Details of the parallel communication are provided in Appendix 7.

Far-field X-ray energy deposition

The X-ray energy deposition module is crucial to the calculation accuracy of the blow-off impulse in the Xablation2D code. To validate this module, we compare the energy deposition rate obtained from the Xablation2D code with that obtained from the Monte Carlo code Geant4 [36-39]. In the Xablation2D simulation, parallel soft X-rays with an initial energy flux ${\varphi }_0=418\text{J}/{\text{cm}}^2$ are incident perpendicularly on a 2D planar Al target with a density $\rho =2.738\text{g}/{\text{cm}}^3$, and the pulse duration is 50 ns. The energy spectrum of the X-rays approximates a blackbody spectrum with a radiation temperature $T=1\text{keV}$. The wavelength range of the energy spectrum is λ=0.1 Å to 10 Å discretized into 23 energy groups for the numerical simulation. The mass absorption coefficients μ in the Al material for the X-rays in each energy group are listed in Table 1.

In the Geant4 code, a parallel blackbody spectrum photon beam with a radiation temperature of 1 keV is configured to simulate soft X-rays. The interaction between photons and materials is simulated using the Livermore low-energy electromagnetic physics model, which is effective within the energy range of photons from 250 eV to 1 GeV. By tracking the trajectory of photons within the Al target, it is possible to quantify the amount of the deposition energy in the material. The energy deposition rates calculated using the two codes are shown in Fig. 3, where the vertical axis corresponds to the energy deposition rate, which signifies the lost ratio of the incident energy flux after traversing corresponding distance in the material, and the horizontal axis represents the depth within the Al target, along with the direction of the incident X-rays. Figure 3 also shows the relative error of the two results by an orange dashed line. The largest deviation is observed near the material surface, which then rapidly decreased to $\mathcal{O}$ (5%) along the incident depth. This is because the energy deposition rate is small at the material surface, where even a minor deviation in parameters such as the mass absorption coefficients in the database can result in a noticeable relative error. In addition, the depth required for 50% of the X-ray energy to be deposited in the material is $4.17\mu \text{m}$ for Xablation2D and $3.5\mu \text{m}$ for Geant4, which is consistent with $\mathcal{O}(5\mu \text{m})$ estimated in [19].

Fig. 3Full-screen View

(Color online) Energy deposition rate varying with the material depth along the X-rays incident direction

Shear flow simulation

We use plane shear flow simulations to verify the numerical diffusion from the algorithms of AWM, IWM, and IIWM. In the Xablation2D code, the size of the simulation domain is $L_x\times L_y=0.4\times 0.4{\text{cm}}^2$ and the domain is discretized into $400\times 400$ grids. 10×10 pseudo-particles are allocated to each grid. The entire flow field is divided into two layers. The upper and lower fluids are denoted by A and B, respectively. The interface between the two fluids is located at $y=0.2\text{cm}$. The mass density, velocity in the x direction, and velocity in the y direction of the upper fluid are ${\rho }_A=2.738\text{g}/{\text{cm}}^3$, $u_{Ax}=3\times {10}^5\text{cm/s}$ and $u_{Ay}=0\text{cm/s}$, respectively. The lower fluid has the same mass density as the upper ${\rho }_B={\rho }_A$, and the velocity is equal in magnitude but opposite in the x direction $u_{Bx}=-u_{Ax}$, and the same velocity in the the y direction $u_{By}=u_{Ay}$. The simulation employs a periodic boundary condition in the x direction, and an open boundary condition in the y direction. In addition, the ideal gas EOS is adopted. $p=(\gamma -1)\rho \mathrm{,}$ (32) where $\gamma =1.667$ is the ratio of the specific heat and the initial pressure p is set to 1 GPa. The fluid velocity distribution in the x direction (ux) of shear flow is shown in Fig. 4.

Fig. 4Full-screen View

(Color online) Fluid velocity distribution in the x direction in the simulation by two different algorithms at various time intervals. (a) AWM; (b) IIWM (or IWM)

In Fig. 4, it can be observed that the IIWM (or IWM) has a better suppression for numerical diffusion than that of the AWM. For the AWM in Fig.4 (a), the jumped shear flow velocity becomes smooth, and this smoothness gradually saturates into the upper and lower layers owing to numerical diffusion, but for the IIWM (or IWM) in Fig.4 (b), the jumped velocity can be maintained throughout the entire simulation time because of using the integration of the internal distribution function in the shear velocity direction.

Kelvin-Helmholtz instability simulation

When two contiguous fluids flow with shear velocities, instability can arise at the interface between the fluids if there is a small fluctuation. This phenomenon is known as the Kelvin-Helmholtz instability (KHI) [40, 41]. The KHI is widely observed in nature, including the turbulent mixing of fluids in jet streams, the formation of undulatus, and supernova explosions [42-44]. Here, we perform KHI simulations using the Xablation2D code and the open-radiation magnetohydrodynamic simulation code FLASH4 [45], and compare the two simulation results. The initial parameters, including the size of the simulation domain, number of grids, and number of pseudo-particles for each grid, are the same as those in Sect. 4.2, except for the temperature. The initial material temperature was set as $T=21.6\text{eV}$, which corresponds to a specific internal energy $e=1\times {10}^5\text{J/g}$. The initial pressure determined using the EOS is 182.24 GPa. Initially, a velocity perturbation in the y direction is introduced at the interface of the two fluids as a seed for the KHI in the following form [46] $u_y^1=u_0\sin \left(kx-\frac{\pi }{2}\right)e^{\left(k|y-\frac{Ly}{2}|\right)}\mathrm{,}$ (33) where $u_0=1\times {10}^6\text{cm/s}$ is the initial amplitude and $k=4\pi /L_x{\text{cm}}^{-\text{1}}$ is the initial wave number. To observe the secondary unstable structures of KHI, it is necessary to reduce numerical diffusion; therefore, we adopt the IWM in the Xablation2D simulation. In the FLASH4 simulation, we maintain consistency in the initialization parameters and perturbation but use adaptive mesh refinement (AMR).

Figure 5 illustrates the mass density of the evolution of the lower half fluid. The six subfigures in the first row are the results obtained from the Xablation2D code and those in the second row from the FLASH4 code. Simulation results for the two codes are consistent. It can be observed that the disturbance at the fluid interface gradually increases, and eventually, the two fluids mix and show vortex structures. The IIWM is also employed to simulate the KHI process with the same parameters as in Xablation2D. The simulation results are consistent with those of FLASH4, except for some secondary structures, as shown in Fig. 6. This deviation can be attributed to the relatively higher numerical diffusion in the direction of convective velocity in the IIWM.

Fig. 5Full-screen View

(Color online) Mass density distribution of the lower half fluid at various time intervals in the Xablation2D and FLASH4 simulations, respectively. The first row corresponds to the result from Xablation2D, and the second row corresponds to the result from FLASH4

Fig. 6Full-screen View

(Color online) Mass density distribution of the lower half fluid at various time intervals in the Xablation2D simulation by IIWM

The mixing width η is defined as the difference between the maximum and initial positions of the lower half of the fluid as shown in Fig. 7 (a). The evolution of mixing width is shown in Fig. 7 (b), where the blue and cyan dashed lines correspond to the results obtained from the Xablation2D, and the magenta solid line corresponds to the FLASH4. The dotted lines indicate the relative errors between the two codes, which remain $\mathcal{O}(5\%)$, indicating a certain level of consistency between them. This minor deviation validates the reliability of the Xablation2D code. It is observed that the growth in the mixing width is significantly smaller than that predicted by the classical linear theory of the KHI. This is because the compressibility of the simulation and evolution of the fluid are nonlinear owing to the strong initial perturbation.

Fig. 7Full-screen View

(Color online) Evolution of the mixing width and relative error in the Xablation2D and FLASH4 simulations, respectively. (a) The mixing width definition; (b) The mixing width at various time intervals

Vaporization blow-off impulse

When X-rays irradiate a solid material, a blow-off impulse load is applied to the material surface, as shown in Fig. 8. The energy of X-ray photons is primarily absorbed by the material through photoelectric effects and Compton scattering. The X-ray energy flux decreases exponentially as it penetrates the interior of the material from the surface. Consequently, the specific internal energy of the material at the surface increases rapidly, which leads to a phase transition and adiabatic expansion, and finally generates an ablation layer. If the deposited energy exceeds the sublimation energy of the material, the solid material at the surface changes into a gas, forming a vaporization zone on the front. The vaporized material is violently sprayed into the vacuum to generate a recoil impulse loading on the remaining material, where the recoil impulse is known as the blow-off impulse (BOI).

Fig. 8Full-screen View

Simple illustration of the far-field X-ray ablation. The shaded area represents the blow-off impulse, and the black solid curve is the energy deposition curve

In the BOI simulation case, the Al material is selected as the target material, the IIWM algorithm is employed in the Xablation2D code. The expansion and compression of the material are described by the PUFF EOS (Eq. 10) and Güneisen EOS (Eq. 12)as discussed in Sect. 2.3. The physical properties of the Al used in the simulation are listed in Table 2. The size of the simulation domain is $L_x\times L_y=0.8\text{cm}\times 0.8\text{cm}$, which is divided into $400\times 400$ grids. We allocated $100\times 100$ pseudo-particles to each grid near the surface of the target material and 5 ×5 pseudo-particles in other regions. Two geometric configurations of the target were simulated: a semi-infinite slab and a cylinder.

4.4.1

Semi-Infinite Slab

The semi-infinite slab is located in the region where x > 0.4 cm, whereas the region where $x⩽0.4\text{cm}$ is filled with low-density gas to maintain the numerical stability of the code. A periodic boundary condition is applied in the y direction, and an open boundary condition is applied in the x direction. Parallel X-rays have a radiation temperature T = 1 keV and are incident on the slab in the normal direction. The energy spectrum was discretized into 23 energy groups, as shown in Table 1. The pulse duration and initial energy flux are 50 ns and 418 J/cm², respectively. Figure 9 illustrates the evolution of the mass density and pressure. The material surface is vaporized by X-ray irradiation and then expands outward into the surroundings. Material expansion generates a blow-off impulse (BOI) and produces a thermal shock wave that propagates inward. The density accumulation in the left region is a consequence of the presence of low-density gas during initialization.

Fig. 9Full-screen View

(Color online) Mass density and pressure distribution in the semi-infinite slab at various time intervals. (a) Mass density varying with the simulation time; (b) The pressure varying with the simulation time

For the semi-infinite slab configuration, the impulse distribution was uniform in the y direction. It is possible to obtain the evolution of the one-dimensional average BOI by integrating it along the y direction and dividing it by the length Ly. This result is indicated by the dashed blue lines in Figs. 10. To verify its correctness, we develop a one-dimensional Lagrangian code, referred to in Ref. [47] to calculate the BOI under the same parameters. The result is shown as a solid magenta line, which is consistent with the Xablation2D code. In addition, we compare the simulation results with the three BOI models, shown as dotted lines in Fig. 10. Their analytical formulas are as follows [9]: $\begin{array}{c}\text{BBAY: }I\text{=}\alpha \sqrt{2}{\left\{{\displaystyle {\int }_0^{x_0}}\left[e\left(x\right)-e_{\text{s}}\right]{\rho }_{0}x\text{ d}x\right\}}^{1/2}\mathrm{,}\\ \text{Whitener: }I\text{=}\sqrt{2}{\displaystyle {\int }_0^{x_0}}{\left[e\left(x\right)-e_{\text{s}}\right]}^{1/2}{\rho }_0\text{ d}x\mathrm{,}\\ \text{MBBAY: }I\text{=}\alpha \sqrt{2}{\left\{{\displaystyle {\int }_0^{x_0}}e\left(x\right)-e_{\text{m}}\left[1+\ln \frac{e\left(x\right)}{e_{\text{m}}}\right]{\rho }_0^{2}x\text{ d}x\right\}}^{1/2}\mathrm{,}\end{array}$ (34) where $e(x)$ is the energy deposition profile, e_m is the melting energy, α is a correction parameter ranging from 1 to $\sqrt{2}$, and x₀ is the thickness of the sublimation layer determined by setting $e(x)=e_{\text{s}}$. In the simulation, we set $e_{\text{m}}=3.975\text{kJ/g}$ and $\alpha =1.1$, respectively. The BOI models estimated only the instantaneous deposition energy for irradiation, resulting in the horizontal dotted lines in Fig. 10. The BOIs calculated using the Xablation2D and 1D Lagrangian codes are shown as blue dashed and magenta solid lines in Fig. 10, has a growth time of $\mathcal{O}(0.1\mu \text{s})$. It demonstrates a tendency to stabilize, which corresponds to the characteristic time in which the material completely sublimates at the surface. The results obtained from the BBAY model and MBBAY model align well with the stable values obtained from the numerical simulations.

Fig. 10Full-screen View

(Color online) Blow-off impulse varying with the simulation time

Furthermore, we compare the simulation results of the Xablation2D code with the published experimental results to verify the reliability of the code. The Ref. [48] provides three measurements of the BOI produced by X-rays irradiating a flat Al material. The X-ray parameters and measured impulse values in experiments are presented in Table 3. We simulated these three experiments using Xablation2D with experimental parameters. The BOIs used in the simulations are shown in Fig. 11(a). The stable values of these BOIs are 97.84, 121.82, and 125.36 Pa⋅s, respectively. Figure 11(b) shows the comparison between the simulation results and the measured BOIs. It can be observed that there are two consistent BOIs, and the relative errors between the simulation and experimental results are negligible for No. 01154 and No. 01170. Although the relative error in case No. 01171 is larger, it remains within $\mathcal{O}(20\%)$.

Table 3Full-screen View

X-ray parameters and measuring impulse in experiments

Number	01154	01170	01171
Radiation temperature T (keV)	0.21	0.227	0.211
Pulse duration ${\tau }_0(\text{ns})$	53	36	44
Initial flux ${\text{Φ}}_0({\text{J/cm}}^2)$	163	181	192
Measuring impulse $I(\text{Pa}\cdot \text{s})$	99.1	118.5	162.2

Show more

Fig. 11Full-screen View

(Color online) Comparison between the Xablation2D simulations and experiments. (a) Blow-off impulse varying with the simulation time; (b) Simulation results vs. the experiment results

4.4.2

Cylinder

In the cylindrical configuration, the target is positioned at the center of the simulation domain. The central coordinates of the cylinder are $x=y=0.4 \text{cm}$, and the radius of the target is 0.2 cm. The remaining areas of the simulation domain are filled with a low-density gas. The physical properties and X-ray parameters used in this configuration are identical to those employed in the semi-infinite slab configuration. The mass density and pressure distributions in the simulation are shown in Fig. 12. X-rays irradiate the left side of the cylinder. The vaporized material is produced on the left surface, ejected violently, and generates a BOI that drives a thermal shock wave moving toward the center of the cylinder.

Fig. 12Full-screen View

(Color online) Mass density and pressure distribution in the cylinder at various time intervals. (a) Mass density varying with the simulation time; (b)The pressure varies with the simulation time

When a parallel X-ray irradiates a curved surface, the density of the deposited energy on the material surface is nonuniform. The BOI is a function of the radial direction of the cylinder. Refs. [11, 10, 49, 21] propose that if the energy deposition density has a cosine profile on the cylinder surface, the variation of the BOI is proportional to the cosine of the polar angle $I=I_0\cdot \cos \theta \mathrm{,}$ (35) where θ denotes the polar angle defined in Fig. 13(a), I₀ represents the BOI at the $\theta {=0}^{\circ }$.

Fig. 13Full-screen View

Blow-off impulse varying with the polar angle at various time intervals. (a) Definition of the polar angle; (b) t=0.04 μs; (c) t=0.10 μs; (d) t=0.16 μs

Figure 13(b), (c), and (d) show the distribution of the BOI via the polar angle at different time intervals. The BOI profiles satisfies the cosine law; however, small deviations also exist. The main reason for this is that the thickness of the X-ray energy deposition layer cannot be completely ignored. The deviation between the BOI curve and the cosine function increases as the polar angle increases from $0^{\circ }$ to ${90}^{\circ }$, particularly at $\theta {=90}^{\circ }$, where the specific internal energy of the material is not zero because of the penetration of X-rays into the material; thus, the BOI is also not zero. The simulation results confirm this phenomenon.