Characteristic method for isentropic compression simulation

Quasi-isentropic compression^1–6 is one method to achieve a high-density, low-temperature state in the laboratory; this is an important research area in the study of astrophysics and the fundamental physics of condensed matter.^7–10 Experimental design and data processing are key components of this experiment. Methods such as hydrodynamic code simulation, Lagrangian analysis,^11–14 and back integration¹⁵ are widely used and have played important roles in this area of study, but there are some deficiencies. The newly-developed iterative characteristic method^16,17 has been the high-performance method for experimental design and data processing method in recent years.

A well-developed hydrodynamic code can give an excellent forward calculation for an isentropic compression experiment,¹⁸ but it can only be operated in the forward direction. It requires iterative calculations to obtain an optimal loading pressure profile. The Lagrangian analysis method is the earliest and simplest data processing method. In situ particle velocities on two or more interfaces are required for this method. To measure the in situ particle velocity, there should be at least two different layers in the target: the sample and the observation window. Researchers always hope that the impedances of these two layers strictly match; however, generally speaking, although great efforts have been made, the impedances cannot match exactly, especially when the pressure is raised. In addition, it is impossible to always find a suitable impedance matching window for every sample, not to mention the difficulties in target fabrication. To solve this problem, D. Hayes and S. D. Rothman have developed the back integration and the iterative characteristic method, respectively, which do not require impedance matching. These two methods are used in the backward direction in a single layer when the rear surface velocity of this layer has been measured. Although more precise stress-density curves can be obtained, there are still some deficiencies in these two methods. A time differential exists in hyperbolic equations with the back integration method, and there are usually some sources of noise in the experimental data. When the time-cell size is smaller than the time-scale of the noise, the time differential may be a huge value, which is a deviation from the physical nature. Furthermore, this error can be carried to the next space grid and may be amplified. Numerical smoothing has been used in this method, but some useful physical information may be lost. The iterative characteristic method avoids this limitation, and it has played important role in data processing even experimental design. If the interaction between characteristics and interface can be properly described, this method can be operated in two directions. Subsequently, experimental design, experimental simulation, and data processing all can be performed by this method.

In this study, we developed such a characteristic method using a Murnaghan-form isentropic equation.^19–22 Herein, first, the characteristic reverberation and transmission on an interface and characteristic intersection are analytically described. Then, an example description of the calculation of a two-layer flow field is given. Finally, example uses of this method are presented, and features of this method are discussed.

To describe this method clearly, we introduce characteristics and Riemann invariants first. When a flow field is regarded as inviscid and isentropic, the equations of mass conservation and momentum conservation can be expressed in a simple form^23–25 

Where x, t, ρ₀, p, u, and c_L are position, time, normal density, pressure, particle velocity, and sound speed, respectively; ± are the propagating directions, representing forward and backward propagation, respectively. The solutions of the equation satisfy du ± dp/(ρ₀c_L) = 0 along dx/dt = ±c_L. The integrals of du ± dp/(ρ₀c_L) are called Riemann invariants

and the routes dx/dt = ±c_L are defined characteristics.^23–25 Regardless of whether the characteristics are bending or straight, the Riemann invariants remain unchanged along the characteristics. Obviously, if the Riemann invariant along every characteristic is known or can be derived, the whole flow field can be determined.

Apparently, if the pressure along an isentropic curve is a single-valued function of density, the Riemann invariants can be simplified to a function of u and c. We use such an equation, which is the isentrope of the simplified equation of state for condensed matter;^19–22 this simplified equation of state has a form similar to the Mie-Grüneisen equation and is expressed as

where e, c₀, and ρ are internal energy, normal sound velocity, and density, respectively. When there is no entropy increase, dS = 0, i.e., de/dρ = p/ρ². The integration form of this equation is

Here, β = (γ − 1)K = β(S). β(S₀) can be derived by setting p = 0. Then, the isentropic line can be expressed by

This equation is similar in form to the Murnaghan equation. With a small entropy increase, this equation has been proven to remain valid, but the coefficient γ₀ should be replaced by a common γ. Then, the isentropic equation is expressed as

The Riemann invariant expressions can be deduced using Eq. (2) and Eq. (6)

The next task is to find the Riemann invariants along every characteristic. The Riemann invariants on the loading surface can be easily calculated from the loading source and primary state of the sample. However, the invariants generated from the interface are more difficult to derive because there are two kinds of samples with different state equations. Additionally, the relationships between regions separated by intersecting characteristics need to be studied. We study the intersection between interface and characteristics first.

To describe this clearly, we study the parameter relationships of a common two-layer target, composed of layer A and layer B, as shown in Fig. 1(a). Layer A is a sample and layer B is an observation window. We choose the Lagragian coordinate system and we are only concerned with three representative characteristics: a forward characteristic CD, its reflection DF, and transmission DE. The whole (x, t) plane is divided into five regions, i.e., A1, A2, A3, B1, and B2, by the three characteristics. A1 and B1 are wave front regions, which are known. A2 also can be described when the loading pressure is known. The regions A3 and B2 may be deduced when these three regions are known.

We can obtain the relationships between these regions along or across characteristics. A1 and A2 are related by a backward characteristic, which is not presented in Fig. 1. Along this characteristic, the Riemann invariant in regions A1 and A2 is

and the one between regions A2 and A3 is

The invariant between regions B1 and B2 is

where u_A1, u_A2, u_A3, u_B1, u_B2, c_A1, c_A2, c_A3, c_B1, and c_B2 are particle velocities and sound velocities in those five regions, respectively; c_A0 and c_B0 are the normal sound velocities of the two layers. In addition, the particle velocities and pressures along the interface are continuous, so the parameters between regions A1 and B1 satisfy

Moreover, the parameters between regions A3 and B2 satisfy

where p_A1, p_A3, p_B1, and p_B2 are the pressures in regions A1, A3, B1, and B2 respectively. The sound velocity function in region A3 can be derived using Eq. (8) to Eq. (14).

Then, the particle velocity in region A3 can be obtained using the Riemann invariant expression

The sound velocity and particle velocity in region B2 can be derived as

The method to calculate the parameters along the interface has been given above. We study the relationships between regions separated by intersecting characteristics next. Four representative characteristics and four split regions, i.e., A4, A5, A6, and A7, are shown in Fig. 1(b). The parameter relationships between regions satisfy the following equations

We can simplify the four expressions above into two compact ones

Obviously, if the parameters of any three regions are known, the parameters of the fourth can be determined.

We now provide an example of the calculation of the whole flow field of a representative two-layer target. The flow field is shown in Fig. 2, with the loading surface at x = 0 and interface at x = h. The whole flow field is divided into 13 regions. When the loading pressure on the loading surface is known, we can calculate every region. Here, region 1 and region 10 are wave-front regions, and region 2 is determined by loading pressure history; therefore, these three regions are known. Region 4 and region 11 can be first calculated using Eq. (15)–(18). Then, region 5 can be determined using Eq. (23)-Eq. (24), and the other adjacent regions are calculated sequentially. When the interface velocities are experimentally measured, we also can calculate the whole flow field from the interface until we obtain the loading pressure.

We performed checks on this method by comparing the calculated results with reported experiments. Two representative target configurations are selected. One is a multilayer target without impedance matching, and the other is a single-layer target without an observation window.

The application to the multilayer target is introduced first. D. B. Reisman performed a ramp compression experiment using the Z accelerator on stainless steel.²⁶ The ramp loading pressure duration was approximately 200 ns with a 16-GPa peak value, while an 800-μm steel sample was placed on a copper driver panel, and a LiF window was used as an observation window. The steel/LiF interface velocity was measured using a Visar. We calculated different particle velocities and pressure profiles (shown in Fig. 3) using this method by setting the loading pressure boundary condition equal to the curve given in Fig. 3 of Ref. 27, reproduced in Fig. 3(b). The density of stainless steel was chosen to be²⁶ ρ₀ = 7.90 g · cm⁻³. The other two coefficients in the Murnaghan-form isentropic equation for stainless steel were c₀ = 4.71 μm/ns and γ = 5.02, which were obtained by searching the minimum rms difference between calculational velocity with the experimental one on steel/LiF interface. These three parameters for Cu were chosen to be¹⁹ ρ₀ = 8.93 g · cm⁻³, c₀ = 3.94 μm/ns, and γ = 4.956; for LiF, they were¹⁹ ρ₀ = 2.64 g · cm⁻³, c₀ = 5.15 μm/ns, and γ = 4.4. Apart from the part below 0.05 μm/ns, most information regarding the calculated steel/LiF interface velocity agrees well with the experimental results. The discrepancy here is discussed in section IV. Other calculated results, such as in situ velocities, interface velocities, and interface pressure, can be determined by this method, which are important to experimental design and have been given in Fig. 3. It is thus clear that this method can be successfully used in multilayer impedance mismatching target simulation.

Another application involves calculating J. P. Davis's experiments,²⁷ which were also performed on the Z accelerator. Four free surface velocities of aluminum samples with thicknesses of 920 μm, 1219 μm, 1518 μm, were measured in this experiment. We calculated free surface velocities, loading surface particle velocities, and characteristics using this method by setting the loading pressure boundary condition equal to the curve given in Fig. 2 of Ref. 28, reproduced in Fig. 4(b). All results are shown in Fig. 4. The density of aluminum here was chosen to be²⁷ 2.703g · cm⁻³, and the other two coefficients in the isentropic equation were obtained by fitting Eq. (6) to the experimental stress-density line of aluminum given in Fig. 5 of Ref. 28. c₀ = 5.828 μm/ns and γ = 3.10. The comparison of calculated free surface velocities and experimental ones is presented in Fig. 4. Most surface velocity information can be determined from the calculation, and the calculated results match the experimental ones better for thinner targets, especially for the increasing part. The discrepancies in the increasing part for the thickest target are also discussed in section IV.

In addition to the results given by J. P. Davis, more information about the flow field can be given by this method. Part of this is presented here. Generally speaking, the free surface velocities, even those on the loading surface, have distinct changes at the times when they are reached by the reverberations. We present the calculated results for the loading surface velocities in Fig 4. The velocities vary with target depth, which is useful information in a multilayer target design.

The velocity details presented in Fig. 4 can be reasonably explained by the characteristics of the flow field, which are shown in Fig. 5. We give some example explanations here. The reverberation reached the loading surface earlier in a thinner target, which can increase the forward Riemann invariants. Therefore, the forward Riemann invariants and backward Riemann invariants of the 920-μm target are larger than those of the 1821-μm target. Larger free surface velocities were obtained with a thinner target. Furthermore, the 1821-μm target shows slowly dropping velocity prior to reverberation(reverberation time 485 ns, seen in Fig. 4(a)), while the loading pressure decreased quickly after the pressure peak. The reason for this can be seen using Fig. 5(c). When leaving the loading surface, the characteristics prior to pressure peak converge, but the ones after pressure peak diverge, which can induce stiff rising edge and slow falling edge on free surface velocity.

General information concerning a ramp compression flow field has been well calculated, but some details cannot be described perfectly by this method. For example, there are discrepancies in the part below 0.05 μm/ns for the steel/LiF interface velocity in Fig. 3(a), and in the increasing part of the thickest sample in Fig. 4(b). The are two reasons for these discrepancies. One is that, Lagrangian sound speed in a solid has relation with strength(shear stress), density(strain), longitudinal stress, even phase transformation, especially remarkable in a low pressure loading experiment. So there are humps or hollows on stress-density curves, which can not be completely reflected by the Murnaghan-form isentropic equation for Lagrangian sound speed obtained by this equation is a single-valued function of density. The shear stress of solid material is usually below GPa, so when the loading pressure rises to several MBar, the strength of solid materials can be neglected, and the solid material can be regard as fluid. Then, the discrepancies caused by material strength can be neglected. The other is that, the coefficient γ has been used as a constant in this equation, which is completely acceptable only when there is no entropy increase. But there is some entropy increase in a ramp compression experiment. So errors will be brought when we use this coefficient as a constant. Therefore, the validity of this equation is the key to the accuracy of this method, i.e., whether the Murnaghan-form isentropic equation describes the actual strain curve perfectly.

This method has been used successfully to simulate the ramp compression on multilayer targets by application in the forward direction. General information concerning the ramp compression flow field can be easily calculated once the loading pressure history and the target configuration are given, which is an important development relative to the original characteristic method.^16,17 Obviously, this method is acceptable not only for experimental design, but also for simulation. Further more, this method can also be used to calculate the stress-density curve of a ramp compression experiment. If there are two or more measured velocities in D. B. Reisman's experiment(in section III), we can obtain the stress-density curve by iterative use of this method without using the loading pressure profile.

In summary, a characteristic method was presented. First, the characteristic intersection and reverberation were analytically solved using a Murnaghan-form isentropic equation and Riemann invariants. Then, the validity of this method was verified by comparing calculated results with reported experimental results. Finally, the deficiencies and significance of this method were discussed. The primary applications show that this method can be used to simulate or design an isentropic compression experiment, or even to obtain stress-density curve of an isentropic compression experiment.

This work is supported by Nature Science Fund of China (Grant No. 11005097), the Foundation of China Academy of Engineering Physics (Grant No. 2010A0102003).