1 Introduction

The constitutive model for materials is the basis of numerical simulation to explain the macro phenomena and predict the mechanical responses of structures, and always can be expressed through the mathematical expression of the material stress-strain relationship. Up to now, constitutive modeling is still a hotspot in mechanical research of solid materials, including the polymer bonded explosives (PBXs) as for their evident nonlinearity and asymmetry of the tensile and compression curves^[1-3]. The tension compression stress-strain curves of PBXs are always S-shaped curves. In previous studies ^[4-5], the tensioncompression constitutive models of PBXs were processed independently. Actually, accurately describing the tension-compression nonlinear constitutive behavior of PBX is still a major challenge in mechanical research of explosives. On one hand, the tension-compression nonlinear constitutive model of PBX has to exhibit satisfactory mathematical description accuracy. On the other hand, model parameters should have explicit physical significance and their initial values should be determined easily.

Research on the PBX constitutive model can be dated back to 1984. Browning investigated the PBX constitutive model in 1984 ^[6] and optimized this model in 1989 ^[7]. However, this constitutive model is only a one -dimensional model. The statistical microcrack model (SCRAM) is a well-known model used to analyze the constitutive behavior of solid brittle materials in recent years ^[8-9]. This model has been extensively used in rocks and concrete. As a constitutive model based on microcrack damages and statistical theory, SCRAM could reflect the physical process of material internal response under stress to a certain extent, which is conducive to determine the micro-deformation mechanism of inner materials. The ViscoSCRAM model ^[10-12] is the most well-known PBX constitutive model and is still being used widely at present. In this model, not only the advantages of the SCRAM model are inherited, but also the viscoelasticity of PBX material is considered. However, the ViscoSCRAM model is disadvantageous and inconvenient to apply as for that too much numerous parameters need to be determined. Besides, this model could also not provide a satisfactory practical application effect, and its practicability still needs to be improved. Though many other works on constitutive behavior of explosives has been reported ^[13-15] recently, but as reported by Shunk in 2014, "even with all of the exploratory works, the constitutive models for PBXs are still lacking." ^[16] In this study, a type of quasi-static tension-compression nonlinear constitutive model is discussed based on our early work ^[17], and its parameters are determined according to test data. Numerical implementation of this model is accomplished in the finite element software, then the accuracy of this model is evaluated by the comparison between Brazilian disk tests and their simulation results.

2 Experiment 2.1 Materials and Instruments

PBX-901, a typical [1, 3, 5-triamino-2, 4, 6 trinitrobenzene] (TATB) -based PBX composed of TATB explosive crystal (>94% mass ratio) and fluorous rubber (< 6% mass ratio) in China, was used as the research object in this study. Uniaxial tension and compression tests were performed on a Instron5582 material testing machine. Tensile dumbbell samples (Φ15 mm × 65 mm) and compressed cylinder samples (Φ20 mm × 20 mm) met the associated requirements of the GJB772A-1997 standard of China.

2.2 Test Method and Data Processing

The uniaxial tension tests and uniaxial compression tests were conducted at 20, 25, 30, 35, 40, 45 ℃ and 50 ℃ respectively. Three samples were tested in parallel at each temperature. The test method met the GJB772A-1997 standard of China.

Acquired test data were processed in two steps. First, test data including stress and strain (positive for tension and negative for compression) were recorded. The curve from zero stress to the maximum tensile/compressive stress (numerical value regardless of the signs) was viewed as the stress-strain curve of the tension/compression segment. Then, the tension and compression stress-strain curves were drawn in the same coordinate system. Figure 1 shows the whole stress-strain curves of PBX-901 considering tension and compression at 20 ℃. According to the proposed stress-strain curves, the arithmetic mean values of tensile strength (σ_t), compression strength (σ_c), and elastic modulus close to zero strain (E₀) of the three PBX-901 samples under the same experimental environment were calculated. Curves at other temperatures have a similar shape with those shown in Fig. 1, but different numerical values.

The results at different temperatures are listed in Table 1. From the Table 1, we can get that all the mechanical constants (tensile strength, compression strength and Young′s modulus) decrease with increasing the test temperature.

3 Characteristics of the Stress-strain Curve and Constitutive Model 3.1 Characteristics of the Stress-strain Curve

Figure 1 shows in complete sigmoid shape stress-strain curves of PBX-901 which can be divided into three stages mathematically, namely, low stress, middle stress and high stress stages. In the low stress stage, the stress-strain curve is nearly a straight line and the tangent modulus which is defined as stress variation / strain variation is approximately a constant higher than zero. In the middle stress stage, the tangent modulus decreases gradually, and the maximum bearable stress of the sample occurs when decreasing to zero. In the high stress stage, stress decreases gradually while strain increases continuously, that is, the tangent modulus further decreases into negative until the final sudden failure. Compression process has all the three stages while tensile process only includes low stress stage.

3.2 Deduction of the Constitutive Model

A S-shaped function is need to describe the stress-strain curves. Boltzmann sigmoid, Gaussian cumulative, and Lorentz cumulative are three typical S-shaped curve functions often used to describe the transition behavior of a physical parameter. These functions have been extensively applied in medicine, biology, agriculture, physics, and chemistry^[18-19]. Their mathematical expressions are shown in Eqs.(1), (2) and (3) respectively. In these equations, x is an independent variable and y is a dependent variable. The parameters of the three equations have different physical significance. Four coefficients, namely, Y₁, Y₂, x₀, and dx, represent the theoretical lower limit, theoretical upper limit, independent variable midpoint, and slope correlation coefficient of function y, respectively. Slope value at x₀ is (Y₂-Y₁) /4dx. The three functions used different mathematical forms to describe the same S-shape curve (Fig. 2). In this study, the Boltzmann function with simpler form and wider application was used to deduce the constitutive model.

$ y\left( x \right) = \frac{{{Y_1} - {Y_2}}}{{1 + {{\rm{e}}^{\left( {x - {x_0}} \right)/{\rm{d}}x}}}} + {Y_2} $

(1)

$ y\left( x \right) = \frac{{\left( {{Y_2} - {Y_1}} \right).\arctan \left( {\left( {\left( {x - {x_0}} \right)/{\rm{d}}x} \right) + \left( {{\rm{ \mathsf{ π} /2}}} \right)} \right)}}{{\rm{ \mathsf{ π} }}} + {Y_1} $

(2)

$ y\left( x \right) = \frac{{{Y_2} - {Y_1}}}{2}\left( {1 + \frac{2}{{\sqrt {\rm{ \mathsf{ π} }} }}\int_0^{\frac{{x - {x_0}}}{{\sqrt 2 {\rm{d}}x}}} {{{\rm{e}}^{ - {t^2}}}} {\rm{d}}t} \right) + {Y_1} $

(3)

A four-parameter Boltzmann constitutive model (Eq.(4)) could be generally derived from Eq.(1) as shown in Eq.(4), σ and ε are the stress and strain, respectively (positive for tension and negative for compression); σ_tll, σ_tul, β, and α are the theoretical lower stress limit, theoretical upper stress limit, correction factor at zero stress-strain point, and correlation coefficient of the modulus, respectively; and β and α are dimensionless numbers. Undetermined coefficients in this model have explicit physical significance.

$ \sigma = \frac{{{\sigma _{{\rm{tll}}}} - {\sigma _{{\rm{tul}}}}}}{{1 + {{\rm{e}}^{\alpha \varepsilon + \beta }}}} + {\sigma _{{\rm{tul}}}} $

(4)

According to the curve characteristics shown in Fig. 2 and the physical significance of parameters, the initial values of the undetermined coefficients in Eq.(4) could be derived easily through the following methods:

(1) σ_tll is selected as the compression strength σ_c, and the initial value of σ_tul is higher than the tension strength σ_t, which is a virtual value that is beyond reach. σ_t can be taken as -σ_c.

(2) The zero point correction coefficient β is used to ensure that zero strain corresponds to zero stress. Based on Eq. (4), the initial value (β₀) can be expressed as follows:

$ {\beta _0} = \ln \left( { - \frac{{{\sigma _{{\rm{tll}}}}}}{{{\sigma _{{\rm{tul}}}}}}} \right) $

(5)

(3) According to the definition of the modulus, the expression of the elastic modulus (Eq. (6)) could be derived from Eq. (4). Based on Eq. (6), the initial value of the modulus correlation coefficient of zero strain (α₀) can be derived (Eq. (7)) as follows:

$ E = \frac{{{\rm{d}}\sigma }}{{{\rm{d}}\varepsilon }} = - \frac{{\alpha {{\rm{e}}^{\alpha \varepsilon + \beta }}\left( {{\sigma _{{\rm{tll}}}} - {\sigma _{{\rm{tul}}}}} \right)}}{{{{\left( {1 + {{\rm{e}}^{\alpha \varepsilon + \beta }}} \right)}^2}}} $

(6)

$ {\alpha _0} = \frac{{{E_0}{{\left( {1 + {{\rm{e}}^\beta }} \right)}^2}}}{{{{\rm{e}}^\beta }\left( {{\sigma _{{\rm{tll}}}} - {\sigma _{{\rm{tul}}}}} \right)}} $

(7)

where E₀ is the elastic modulus value close to the zero point.

In particular, if the tension and compression stress-strain curves of the material are symmetric, that is, σ_tul is equal to σ_tll in numerical value and has opposite signs (σ_tll = -σ_tul = σ_c and β = 0), then the mathematical expression can be simplified into the two-parameter Boltzmann constitutive model ex-pressed in Eq.(8). σ_c can be tested directly and the undetermined coefficient α can be fitted or calculated by Eq. (9)

$ \sigma = \frac{{2{\sigma _c}}}{{1 + {{\rm{e}}^{\alpha \varepsilon }}}} - {\sigma _c} $

(8)

$ {\alpha _0} = - \frac{{2{E_0}}}{{{\sigma _c}}} $

(9)

3.3 Establishment of the Constitutive Model

As shown in Figs. 3 and 4, test data of PBX-901 at 20 ℃ and fitting results with the four-parameter Boltzmann constitutive model (Eq.(10)) and the two-parameter model (Eq.(11)) were compared. The mean -adjusted R² values between test data and fitting data in Figs. 3 and 4 are 0.99852 and 0.99712, respectively.This result shows both the fitting curves with both models are highly consistent with the test data, and the description accuracy with four-parameter model (Fig. 3) is slightly higher.

$ \sigma = - \frac{{69.38}}{{1 + {{\rm{e}}^{4.669\varepsilon - 0.441}}}} + 42.26 $

(10)

$ \sigma = - \frac{{55.16}}{{1 + {{\rm{e}}^{5.15\varepsilon }}}} + 27.58 $

(11)

According to the analysis presented in Section 3.2, σ_c in the two-parameter Boltzmann constitutive model (Eq.(8)) can be obtained directly from the test (Table 1). According to Eq.(9), α₀ is 534.45 at 20 ℃ using true strain expression, and it is 5.345 when expressed by percentage strain. This result is approximately 3.625% higher than the fitting one (5.518). Actually, the parameters of the two-parameter Boltzmann constitutive model were not fitted, but the analytic results (σ_c = -27.58 MPa and α₀ = 5.345) were adopted, which can also provide satisfactory description accuracy (Fig. 5). When calculating the analytical results of the parameters, σ_c and E₀ are needed, whereas the tensile curve is unnecessary. In other words, the two-parameter Boltzmann constitutive model of PBX-901 considering tension and compression can be established even if no tension test and data fitting have been conducted.

The simulation results of PBX-901 (parameter used fitting values and analytic values) were compared with the test ones (Figs. 6 and 7) at 50 ℃ to verify the applicability of the two-parameter constitutive model at different temperatures. The fitting values of σ_c and α were-17.58 MPa and 5.566, respectively, whereas the analytic values were-17.58 MPa and 5.552, respectively. The difference of α was only approximately 0.252%.

According to the previously presented analysis, the two-parameter Boltzmann constitutive model could accurately describe the compression constitutive behavior of PBX-901. However, this model can only describe constitutive behavior when stress and strain increase (numerical value) simultaneously, but could not describe the phenomenon when stress decreases and strain increases after the maximum compression load. The parameter values of the two parameter Boltzmann constitutive model within 20-50 ℃ are shown in Table 2.

4 Numerical Simulation and Application of the Model 4.1 Secondary Development of the Two-parameter Boltzmann Constitutive Model

The secondary development of the two-parameter Boltzmann constitutive model is based on user programmable features (UPFs) and usermat subroutine of the ANSYS software. Through UPFs, user-defined material, element and failure criterion (for composites) can be realized through programming corresponding subroutine using FORTRAN language. There are three steps to finish a material subroutine. First, the equivalent strain from stress is calculated according to the constitutive model (Eq.(8)). Secondly, the stiffness and Jacobian matrices are calculated according to the equivalent strain. Finally, the secant modulus is used instead of the elasticity modulus in the stiffness and Jacobian matrices to calculate and update current stress. These steps will be repeated in numerical simulation until the end of computation. Specially, the tangent slope is calculated to replace the secant modulus as for that it′s difficult to calculate when strain is close to zero. This program could be applied in the ANSYS software after being compiled.

4.2 Numerical Simulation of the Brazilian Disk Test

To verify the description accuracy of the two parameter Boltzmann constitutive model in practical engineering, the Brazilian disk test was performed and simulated. The Brazilian disk test is extensively used in investigating the mechanical performances of explosive materials. The test in this study was conducted at 20 ℃ by using three Φ25 mm × 10 mm compressed samples. Sample size and test results are presented in Table 3. The test results show that the average size is Φ24.903 mm × 9.997 mm. Moreover, the average values of critical load and load displacement are almost 1856.441 N and 0.444 mm respectively.

According to the symmetric relationship between stress and boundary conditions in the Brazilian disk test, a simulating model with 1/2 thickness and 1 / 4 perimeter was established. The mean sample size and model parameters at 20 ℃ were used as the input. The SOLID185 element, which is a 8-node finite element, was applied during modeling. During numerical simulation, the tested loading displacement was used as input. Thus, the description accuracy of the two-parameter Boltzmann constitutive model in practical engineering was reflected by the difference between calculated node counterforce of the pressing plate and the tested critical load. According to the symmetric relationship of the model, the plate displacement input in numerical simulation was 0.222 mm. The displacement field, node counterforce, and the first major principal stress output at this moment are shown in Fig. 8. The total node counterforce of the pressing plate is four times that of the model output (440.39 N) and reaches 1761.56 N while the test critical load is 1856.441 N. The relative error is only approximately 5.11% which reflects that the established two-parameter Boltzmann constitutive model has satisfactory description accuracy.

5 Conclusions

A type of quasi-static tension-compression nonlinear constitutive model was deduced by using Boltzmann function in the present work, whose undetermined parameters show explicit physical significance and can be easily determined. Meanwhile, the two-parameter Boltzmann constitutive model was further deduced based on the symmetry of curves, and there are two ways to determine the initial parameter values. One is by fitting with test result and the other just using the compression strength and the initial modulus.

According to the test data, the four-parameter and two-parameter Boltzmann constitutive models of PBX-901 were established under different temperatures. Through analysis, both models showed satisfactory description accuracy. Besides, the two-parameter constitutive model shows a small difference between analytic and fitting parameter values, which indicates that just two physical constants, compression strength and modulus value, from conventional compression test are adequate to establish the high-accuracy tension-compression nonlinear constitutive model, and the tension test is no longer necessary. What′s more, based on the UPFs platform of ANSYS software, a material subroutine of the two-parameter model was developed and the numerical simulation of the Brazilian disk test was performed, which used loading displacement as the input and critical load as the output. The simulation result shows only 5.11% relative error compared with the test results. This small error confirmed the high description accuracy of the two-parameter model in practical applications.