Numerical Simulation of Thermal Shock Waves Induced by Pulsed X-ray in C/PF Materials

Based on the orthotropic elastoplastic constitutive model of carbon-phenolic (C/PF) material under three-dimensional strain, a self-developed explicit finite element program TSHOCK3D is developed to simulate the multi-physical coupling process of thermo-mechanical response of C/PF target irradiated by pulsed X-rays. Two typical irradiation conditions are simulated, namely blackbody temperatures 1keV soft X-ray and 3keV hard X-ray with the same energy flux of 300J/cm2 and pulse width of 100ns. The C/PF target presents two different thermodynamic patterns. Compared with 3keV hard X-ray, 1keV soft X-ray has a shallower energy deposition depth and a higher peak pressure of thermal shock waves. Furthermore, the irradiated surface under 1keV X-ray exhibits the phenomenon of vaporization while 3keV does not. Following the compression wave is a tensile rarefaction wave, the peak pressure of which is higher in the case of the 1keV. In addition to the phase transition and shock wave propagation in the direction of irradiation, the tensile fracture phenomenon appears in the side area of the target plate, and the effect is more significant under the irradiation of 3keV.


Introduction
Carbon fiber reinforced polymer composites (CFRP) represented by carbon phenolic (C/PF), as a new type of composite material, has excellent characteristics such as low density, high strength, well designability, etc. Since it came out, it has been widely and successfully used in the field of aerospace and military industry, e.g., the manufacture of missile shell and space shuttle skin [1,2]. Spacecraft are often under extreme external environmental conditions in outerspace and may be subjected to extreme dynamic loads such as high-energy X-ray irradiation [3]. When materials are irradiated by high-energy pulsed X-rays, the photon energy is transformed into internal energy of the material through a series of complex physical and chemical processes. If the energy is high enough, it can cause the surface materials to vaporize and recoil. At the same time, due to the rapid decline of energy deposition from the surface to the inside, it will form a large specific internal energy gradient and thus pressure gradient in the materials. Eventually, these two physical processes jointly lead to the generation of shock waves in the materials, namely, X-ray thermal shock wave [4,5]. For this multi-physical coupling problem, it is difficult to carry out high-energy X-ray irradiation experiments because of the experimental cost. Alternatively, simulated loading experiments such as explosive loading [8], electron beam loading [6,7] were carried out with lower expenditure. However, the equivalence experimental error cannot be completely eliminated due to different physical mechanisms. Numerical simulation has unique advantages in expenditure, modeling flexibility, data acquisition, and parameter optimization. With the improvement of computer performance, it has become more and more important in this research.
However, in the current numerical simulation work, irradiated materials are usually limited to metallic isotropic materials [9], and the processing of the anisotropic constitutive model mainly involves experimental fitting of the parameters of the traditional equation of state [10,11]. Among the existing numerical simulation studies considering anisotropic correction, examples under one-dimensional strain and two-dimensional plane strain are mainly used [12,13]. This paper presents an orthotropic constitutive model under three-dimensional strain conditions, considering the plastic yield criterion and fracture damage criterion. Based on the model, an explicit finite element method program TSHOCK3D is developed. Taking a rectangular target of C/PF material irradiated by pulsed X-ray as an example, the propagation law and damage characteristics of the three-dimensional thermal shock wave in the target plate under 1keV and 3keV irradiation conditions are compared and analysed.

Orthotropic constitutive model of C/PF materials
2.1. Orthotropic constitutive model of C/PF materials C/PF material was prepared by carbon fiber woven into a carbon cloth, then layered and injected into phenolic resin for heating curing process. In this kind of material, phenolic resin is only used as an adhesive, and the carbon fiber carries most of the load. Therefore, the mechanical properties of the C/PF material are closely related to the properties and distribution direction of the carbon fibers. Due to the economic efficiency of industrial production, carbon fibers are usually orthogonal weaved, thus C/PF materials usually have an orthogonal anisotropy in mechanical properties.
For orthotropic materials, the constitutive model is only established in the principal axis coordinate system of the material. Under three-dimensional strain conditions, the constitutive relation of elastic segment can be described by Hook's law: 11 12 [C] is the stiffness matrix, determined by Young's modulus, Poisson's ratio and shear modulus in the corresponding direction. The stress and strain tensor are decomposed into spheres and deviations, and the pressure in the elastic segment of orthotropic materials can be expressed as: in this paper.  (2) is only valid in the elastic section at low pressure, it is necessary to introduce the equation of state to calculate the nonlinear relationship between the pressure and specific volume of the material under high pressure. For the X-ray irradiation problem, the material may have thermal expansion or even vaporization in the surface area of the target plate, while the inside of the target plate will be subjected to compression loading by shock waves. Therefore, PUFF equation of state [13] and Grüneisen equation of state are adopted to describe the pressure under different states, respectively:

Since the linear relationship of Equation
where p is pressure, p H is Hugoniot pressure, e H is Hugoniot specific energy, v is specific volumn，  is density, 0  is initial density, e is specific internal energy, Γ is Grüneisen coeffcient,  is specific heat ratio of gas, Under the condition of small strain, the following equation is approximately satisfied： The Taylor expansion of Equation (3) is carried out with respect to  , and the higher order terms of more than three times are ignored. Then the polynomial form equation of state with respect to  is obtained as: By taking the full derivative of Equation (5), we can get the pressure dp in the incremental form at a single time step: Equation (5) and Equation (6) are only applicable to isotropic materials, which are essentially equivalent in obtaining the Murnaghan form equation of state for materials. A 1 =B 1 is equivalent to the bulk modulus of the material. However, Equation (6) does not reflect the anisotropy of the material. At elastic segment, the linear relationship between pressure increment and specific volume increment determines that the pressure is mainly affected by the first order term of specific volume. At the same time, the anisotropic pressure is also affected by the partial strain in the material. Therefore, we modified the Equation (6) So far, a modified equation of state for orthotropic materials is obtained in this paper. The pressure calculation of this equation takes into account the nonlinearity of pressure variation with specific volume and the anisotropic properties of materials.

Tsai-Hill yield criterion and plastic strain calculation
For the CFRP materials, based on the experimental measurement work of a certain type of C/PF by Jiang [12,14], Tsai-Hill yield criterion is adopted, and its basic form under the three-dimensional strain condition is:  are defined as 11   22  33  2  2  2  2  2  2  2  2  2  11  22  33  22  33  11  33  11  22   1  1  1  1  1  1  1 13 and Y 23 are the uniaxial tensile (or compression) strength in 1, 2 and 3 directions and the shear strength in 12, 13 and 23 directions, respectively.

Material failure criterion
In this paper, only the tensile failure caused by rarefaction wave is considered, and the failure criterion is the maximum tensile stress criterion: is tensile failure strength in three directions. It should be noted that, for isotropic materials, the maximum tensile stress criterion is to solve the stress tensor after the similarity transformation of the tensor matrix, then to solve the principal stress and determine the failure. However, for orthotropic materials, the failure determination is based on the stress in the three main axes. The three main axes of the material are the warp direction of carbon fiber, the weft direction of carbon fiber and the lamina stack direction of carbon cloth, respectively. Therefore, after obtaining the stress tensor in the principal axis coordinate system, the diagonal element in the stress tensor can be directly determined.

Results and discussions
The rectangular target of C/PF material is modelled as the irradiated object, the thickness of the target plate is 0.5cm, the width and height are 1.0cm, and the FE model is meshed with the hexahedron element with a uniform size of 0.0125cm×0.025cm×0.025cm in x, y, z directions. There are 33 C/PF material parameters in total. See article [15] for details. X-ray energy spectrum is described by black body spectrum, and its wavelength discrete method and material absorption coefficient calculation method refer to [13].
Pulsed X-ray radiate the square target surface parallel to the x-axis direction. This paper mainly studies the influence of blackbody temperature on the thermal shock wave formation mode and the fracture damage mode of CP/F target plate. Therefore, the incident energy flux 0  (300J/cm 2 ) and The X-ray pulse is assumed as a square wave. The 3D X-ray thermo-machanical response simulation program TSHOCK3D, developed by our research group, is used for calculation, and the physical simulation duration is 1.2 μs .
Generally, the thermal shock waves generated by pulsed X-ray irradiation is jointly formed by the vaporization recoil effect of the material surface and the thermal deformation effect of the energy deposition area. The former is formed by the impact loading of vaporized substances on the remaining solid materials, and it will generate a compression wave propagating along the irradiation direction. The latter is due to the exponentially decreasing distribution of energy deposition in the depth of the irradiated surface. According to the Grüneisen EOS, the pressure in a solid material under constant volume conditions is proportional to the specific internal energy. Under the effect of material inertia, pressure gradients will form compression waves that propagate to the surface and inside of the material. Furthermore, the compression waves generated by the two effects reflect on the front and side free surface and form a superimposed rarefaction wave propagating inward. Considering non-transient loading of pulsed X-ray, the superposition of waves generated at different times cause the waveform to be more complicated. Figure 1 and Figure 2 respectively show the two-dimensional and three-dimensional contour of pressure distributions in the target plate under 1keV and 3keV X-ray radiations at different times. The two-dimensional pressure contour is intercepted on the half-height plane (z=0.5cm) of the target.
As shown in Figure 1, the surface material of the target plate is vaporized, and the vaporized material with high temperature and pressure expands freely and violently. A compression shock wave propagates inward from the surface of the target plate under 1keV soft X-ray irradiation. Following compression wave, there is a superimposed rarefaction wave, which makes the elements into a tensile failure.
As shown in Figure 2, the surface material fails to vaporize. The peak value of shock wave is lower and the FWHM is broader than that in 1keV case. Meanwhile, there have been much more tensile-failed elements in the propagation path of the superimposed rarefaction wave.    The pressure distribution at axis line (y=0.5cm, z=0.5cm) at different times are presented in Figure  3. It can be seen that, the pressure waves of 1keV and 3keV irradiation conditions are approximately triangular waveforms. As the waves propagate, the peak pressure is attenuated and the waveforms are more extended.
For 3keV X-rays, its stronger transmission ability leads to a noticeable energy deposition depth up to the rare free surface at initial moment. The compression waves formed near the front and rear free surfaces are reflected by the free surfaces to form rarefaction waves, and the rarefaction waves near the rear free surfaces are not strong enough to unload the pressure to zero. Refer to Figure 2, there is no spallation phenomenon because peak pressure of the rarefaction wave on the free surface is less than the tensile strength.
For 1keV X-rays, the transmission ability is relatively weak, and most of the energy is absorbed by the material near the front surface, resulting in vaporization and a large energy deposition gradient (corresponding large pressure gradient) near the front free surface. x=0.1375cm x=0.2125cm Figure 4. Pressure-time curves at monitoring points As shown in Figure 4, the pressure-time curves under 1keV and 3keV X-ray irradiation conditions are compared. For each conditions, three monitoring points on the central axis are sampled at the distance of 0.0625cm, 0.1375cm, and 0.1225cm from the front surface. Since the pulse irradiation is set as a square wave with constant power in this article, energy deposition is the only pressure increase mechanism before the pressure wave propagates. The pressure at the three monitoring points increases linearly within the irradiation duration 0