Analysing the impact of renewables on Iberian wholesale electricity market prices using machine learning techniques

1.   Introduction

Polymer bonded explosives (PBXs) are energetic particulate composites ^[1,2,3] widely applied as rocket propellant and main explosive charge in munitions. Different from the common matrix-based heterogeneous composites, three unique characteristics of PBXs should be noted purposely: ($i$) The volume fraction $f$ of the explosive crystals is much higher than that of the binder material ^[3,4,5], and often larger than $90\%$; ($ii$) Young's modulus $E_{b}$ of the binder material is in high contrast with that $E_{p}$ of the explosive crystals, and $E_{p}/E_{b}$ can even reach up to 20,000 at the room temperature and low strain rates ^[6]; ($iii$) They are hazardous media ^[7,8] during consolidation, machining, and transportation processes. These typical features introduce many tough difficulties in theoretical modeling, multi-scale simulation, and experimental tests of their inherent properties ^{[1,2,3,4,6,8,9,10,11]}. By now, evaluating the overall moduli of PBXs efficiently and understanding their fundamental mechanism are still largely open problems to be explored.

Three major strategies applied to investigate the elastic properties of PBXs thus far are the experimental test, molecular dynamics (MD) simulation, and micro-mechanics-based method. Among these, experimental tests are the earliest but still widely adopted one, and various original investigations have been carried out on the macroscopical constitutive relation ^{[1,12,13,14,15]}, failure mechanism ^{[8,10,16,17,18,19]}, influence of explosive and polymeric binding materials ^[3,20], microstructure ^[9,21], interfacial imperfections ^[5], and so forth. It is noted that experimental tests often cover a long time-period, and are conducted with small samples ^[3,10] in view of the hazardous nature of PBXs. Hence, it is pretty time- and money-consuming to investigate the influence of a certain number of factors thoroughly.

MD simulation is another prospective strategy to predict the mechanical properties of PBXs from the molecular level ^[2,22,23,24]. Most current works with this method mainly concentrate on the influence of lattice structures and interactions between different composition elements. However, the effect of explosive crystal size/shape, binder distribution, and other factors can be hardly taken into account by far ascribed to the restriction of the computer hardwares.

Therefore, our attention of this work is drawn primarily on the third strategy, i.e., the micro-mechanics-based method. Within this framework, many creative works have been done. Banerjee and Adams ^[6,25,26,27] developed a recursive cell approach to rapidly evaluate the effective elastic properties of PBX 9501 based on the real-space renormalization group method. Tan and his coworkers ^[5,7,28] established an extended Mori-Tanaka method to compute the effective bulk and shear moduli of PBXs involving both linear and nonlinear interface debonding behavior. Wu and Huang ^[1] first presented a micro-mechanical model by combining a continuum damage model and a rate-dependent cohesive model to characterize the damage of explosive crystals and interface debonding within PBXs, and then proposed a numerical approach framed within the finite element method to study the damage effect induced. Barua and Zhou ^[29,30] established a cohesive finite element framework to quantify the thermomechanical response of PBXs at the microscale. Dai et al.^[4] applied the finite element method to compute the effective moduli of PBXs, and discussed the influence of volume fraction, gradation of explosive crystals, and binder in detail. For the sake of simplicity in theoretical analysis and numerical modeling, most of these works assume that the explosive crystals are in 2D simplified or 3D regular shapes (sphere, hexahedron, etc). However, explosive crystals of practical interest are 3D arbitrarily shaped and randomly distributed according to the experimental observations and optical micrographs (c.f., ^[3,4,10,30]). Furthermore, not only the overall moduli ^[4,5,7,28] but also the physical behaviors ^[5,9,30] of PBXs may be significantly varied by the size and shape of explosive crystals and binders, particularly when interface debonding occurs ^{[4,31,32,33,34]}. In the previous works, a physics-based general imperfect interface model ^[32,33,35] and a physics-based spring-type interface model^[36], which is equivalent to the practical thin interphases for particulate and fiber-reinforced composites containing three or more phases, was developed by combining the Hadamard's relation and rigorous theoretical deduction. This paper aims to further deduce the physics-based imperfect interface model to characterize the thin soft binder interphases of highly-filled two phase PBX media, which will bring much convenience for modeling the explosive particles of arbitrary shapes and developing innovative computational methodology to predict the performance of PBXs, and to design and analytically solve the overall moduli of a simplified benchmark model containing the physics-based imperfect interface model derived.

This paper is organized as follows. In the next section, a spring-type interface model with apparent physical background is derived based on the general interfacial relations in ^[32,33] to characterize the thin binder interphases inside PBXs. Afterward, the averaged together with the effective quantities over a representative volume element (RVE) of PBXs are reformulated based on their primitive definitions in micro-mechanics. Section 4 is then dedicated to designing a simplified PBX model and deriving its effective isotropic moduli as well as the upper & lower bounds of the moduli. Some investigations are conducted on the size dependent effective moduli of PBXs in Section 5, and a few conclusions are drawn eventually.

2.   Reformulation of the binder interphases in PBXs as zero-thickness imperfect interfaces

As mentioned in the last section, the thickness of binder interphases in PBXs is pretty thin and irregularly distributed ^{[3,8,10,18,19,20]}, which introduces much difficulty in theoretical analysis and efficient numerical modeling. Hence, the following assumptions are made on the binder interphases so as to developed the theoretical interface model such that (i) the thickness of each binder segment between two adjacent explosive crystals is uniform (see Figure 1); (ii) the binder segments between different explosive crystals may have different thickness; (iii) each interface between the binder interphase and the adjacent explosive crystal is perfect. In what follows, we shall reformulate such binder interphases based on these assumptions as spring-type interfaces of null thickness, which are physically equivalent to the practical binder interphases.

2.1.   Constitutive equations

Let's first consider the 3D PBX microststructure in Figure 1. The whole domain $\Omega$ contains $m$ explosive particles, namely, $\widehat{\Omega }^{\left(r\right) }$ ($r = 1, 2, \cdots, m$), bonded together by the binder interphases $\Omega ^{\left(0\right) }$. Each binder segment is uniform in thickness, while its thickness value allows it to vary from one another (c.f., Figure 1). The interface between a generic explosive particle $\widehat{\Omega }^{(r)}$ and its adjacent binder interphase $\Omega ^{(0)}$ is regarded as being perfect.

In this study, the materials forming $\widehat{\Omega }^{\left(r\right) }$ and $\Omega ^{\left(0\right) }$ are supposed to be individually homogeneous and linearly elastic according to the references ^[20,27,28]. Hence, their constitutive equation can be expressed as

with $\boldsymbol{\sigma }^{\left(\gamma \right) }$ and $\boldsymbol{ \varepsilon }^{\left(\gamma \right) }$ separately denoting the second-order Cauchy stress and strain tensors defined over $\widehat{\Omega }^{\left(\gamma \right) }$ $\left(\gamma = 0, 1, \cdots, m\right)$. $\mathbb{C} ^{\left(\gamma \right) }$ is the fourth-order elastic tensor over $\widehat{ \Omega }^{\left(\gamma \right) }$, and its components satisfy the symmetric property such that

Experimental tests of pure explosive particle and binder materials ^{[5,6,20,27,28]} indicate that both media are isotropic materials at low strain rates; thus, we mainly focus on this particular and important kind of PBXs in this work, i.e.,

where $\kappa_{b}$ and $\mu_{b}$ denote the bulk and shear moduli of the material forming $\Omega ^{\left(0 \right) }$, and $\kappa_{p}$ and $\mu_{p}$ are the bulk and shear moduli of the material constituting $\widehat{\Omega} ^{\left(r\right) }$. They can be further expressed in terms of the Young's modulus $E_{\iota}$ ($\iota = p, b$) and Poisson's ratio $v_{\iota}$ as

Symbols $\mathbb{J}$ and $\mathbb{K}$ are the isotropic operators and satisfy

Above, $\mathbf{I}$ and $\mathbb{I}$, respectively, represent the second- and fourth-order identity tensor.

Geometrical equation of $\widehat{\Omega }^{\left(\gamma \right) }$ fulfills the relation

where $\boldsymbol{u}^{\left(\gamma \right) }$ refers to the displacement vector over $\widehat{\Omega }^{\left(\gamma \right) }$, and $\nabla$ is the gradient operator.

Thus far, all constitutive equations of each phase $\widehat{\Omega } ^{\left(\gamma \right) }$ of the PBX microstructure have been introduced and will be applied to establish the equivalent description of the binder interphase in the next subsection.

2.2.   Reformulation of binder interphases based on the asymptotic analysis

We now consider a practical situation where two generic explosive particles $\widehat{\Omega }^{\left(1\right) }$ and $\widehat{\Omega }^{\left(2\right) }$ are bonded together by a thin binder interphase $\Omega ^{\left(0\right) }$ of thickness $h$ (see Figure 2(a)). Precisely, $\Omega ^{\left(0\right) }$ is assumed to be in a curved configuration for a universal purpose. The middle interface of $\Omega ^{\left(0\right) }$ is donted as $\Gamma _{0}$, and its unit normal vector is specified as $\boldsymbol{n}$. Interfaces $\Gamma _{1}$ between $\Omega ^{\left(0\right) }$ and $\widehat{\Omega }^{\left(1\right) }$, and $\Gamma _{2}$ between $\Omega ^{\left(0\right) }$ and $\widehat{\Omega } ^{\left(2\right) }$, are regarded as being perfect.

Since $\Omega ^{\left(0\right) }$ brings many tough troubles for both theoretical analysis and numerical modeling, in what follows, we shall substitute $\Omega ^{\left(0\right) }$ by an interface $\Gamma _{0}$ of null thickness. Meanwhile, the materials forming $\widehat{\Omega }^{\left(1\right) }$ and $\widehat{\Omega }^{\left(2\right) }$ are extended to $\Gamma _{0}$; see Figure 2(b). The newly generated domains are named as $\Omega ^{\left(1\right) }$ and $\Omega ^{\left(2\right) }$, respectively. In our previous works ^[32,33], the general equivalent interfacial relations across $\Gamma _{0}$ in Figure 2 have been derived by applying the Taylor expansion technique and the continuous condition across $\Gamma _{1}$ and $\Gamma _{2}$ such that

where $\boldsymbol{N} = \boldsymbol{n\otimes n}$ and $\boldsymbol{T} = \boldsymbol{I}-\boldsymbol{N}$. Symbol $[\left[ \boldsymbol{\bullet }\right] ]$ is the jump operator defined by $[\left[ \boldsymbol{\bullet }\right] ] = (\boldsymbol{\bullet })^{(+)}-(\boldsymbol{\bullet })^{(-)}$, and $\langle \bullet \rangle = 0.5\left(\bullet ^{(+)}+\bullet ^{(-)}\right)$ refers to the interfacial average operator. $(\boldsymbol{ \bullet })^{(+)}$ denotes the quantity $(\boldsymbol{\bullet })$ value evaluated at the interface $\Gamma _{0}$ on the side of $\Omega^{\left(2\right) }$ where the normal vector $\boldsymbol{n}$ is directed into; $(\boldsymbol{\bullet })^{(-)}$ represents the quantity $(\boldsymbol{\bullet })$ value also at $\Gamma _{0}$ but from the side of $\Omega^{\left(1\right) }$ where $\boldsymbol{n}$ is originated from. $h$ stands for the thickness of the binder interphase replaced by the energetic interface $\Gamma _{0}$. Normal traction vector $\boldsymbol{t}$ is defined through $\boldsymbol{t} = \boldsymbol{\sigma \cdot n}$, and surface strain $\boldsymbol{ \varepsilon }_{s}$ can be computed by $\boldsymbol{\varepsilon }_{s} = \boldsymbol{T}\cdot \boldsymbol{\varepsilon }\cdot \boldsymbol{T}$. $div_{s}(\mathbf{\bullet })$ is the surface divergence fulfilling $div_{s}(\boldsymbol{\bullet }) = \nabla (\boldsymbol{\bullet }):\boldsymbol{T}$. Parameters $c_{1}\sim c_{5}$ are the material-related coefficients specified by

with $b_{1}^{(r)} = \left(3\kappa ^{(r)}-2\mu ^{(r)}\right) /\left(3\kappa ^{(r)}+4\mu ^{(r)}\right)$ and $b_{2}^{(r)} = 3/\left(3\kappa ^{(r)}+4\mu ^{(r)}\right)$ defined over $\widehat{\Omega }^{\left(r\right) }$ ($r = 1, 2, 0$).

It is noted that the error of the physical fields (i.e., displacement and stress fields) inside the practical three-phase and the foregoing equivalent two-phase models are kept within a fixed bound of $\mathbf{0}(h^{2})$. Interested readers are referred to ^[32,33] for more detailed deductions and analysis.

Next, we invoke the fact that: (ⅰ) All explosive crystals have the same mechanical property, thus, $\kappa ^{(1)} = \kappa ^{(2)} = \kappa_{p}$ and $\mu^{(1)} = \mu^{(2)} = \mu_{p}$. (ⅱ) The stiffness of binder interphases is much lower than that of the explosive crystals ^[6,27], i.e., $E^{(0)}\ll E^{(2)}$. To facilitate the deductions, a small dimensionless parameter $\eta$ is introduced such that $\eta = h/h_{0}\ll 1$, with $h_{0}$ being the reference length of the same order as the geometrical dimensions of the explosive particle in PBXs. Accordingly, the material property of the interphase may be expressed as

where $\hat{\kappa} ^{(0)}$ and $\hat{\mu}^{(0)}$ are the reference binder stiffness parameters of the same order with $\kappa_{p}$ and $\mu_{p}$.

Taking account of the relations

and introducing Eq (2.9) into (2.6) and (2.7) delivers the following compact interfacial relations of the binder interphases:

with

Equations (2.11) and (2.12) indicate that the displacement field presents a jump across $\Gamma _{0}$, but the normal traction field remains continuous within the whole domain. These two interfacial relations are also known as the spring-layer interface model in open literatures (see, e.g., ^[31,36]). It is noted that though the foregoing interfacial relations are derived on the basis of the fact that the thickness of each binder segment between two adjacent explosive crystals is uniform, the thickness of each elemental interface can be randomly defined when the numerical methodology, say eXtended Finite Element Method (XFEM) ^[33,36] and so forth, is employed to investigate the binder effect.

3.   Overall moduli evaluation of PBXs using the micromechanics

Macroscopically, PBXs exhibit linear isotropic elasticity ^[4,6,7,12]. Hence, their constitutive formula is governed explicitly by

where $\kappa ^{\ast }$ and $\mu ^{\ast }$ separately denote the effective bulk and shear modulus of PBXs. Symbols $\mathbf{\Sigma }$ and $\boldsymbol{E}$ are the second-order stress and strain tensors at the macroscopic scale, and $\boldsymbol{(\bullet)}^{\prime }$ refers to the deviator part of the tensor $\boldsymbol{(\bullet)}$ defined as

To evaluate the effective moduli of PBXs in Eq (3.1) by using the micro-mechanics, we take the RVE model $\Omega_{RVE}$ of side length $w$ in Figure 3 as an illustration for the afterward demonstration and deduction. It is noted that the equivalent binder interfaces between energetic crystals can be located at the external surfaces of the PBX RVE, and their influence can be properly considered by taking account of the aforementioned interfacial relations in Eqs (2.11) and (2.12). Precisely, as the RVE model is small enough that it can be viewed as one point from the macroscopic scale and large enough from the microscopic level containing all the micron characteristics of PBXs, the macroscopic stress $\mathbf{\Sigma }$ and strain $\boldsymbol{E}$ tensors can be evaluated equivalently by

with $V$ the volume of $\Omega _{RVE}$, and $\boldsymbol{n} _{R}$ the unit vector normal to the external boundary $\partial\Omega _{RVE}$ of $\Omega _{RVE}$.

In this study, the displacement-controlled boundary condition $\boldsymbol{u} = \boldsymbol{\varepsilon }^{0}\boldsymbol{x}$ is considered and enforced on all external surface $\partial \Omega _{RVE}$ of the RVE model to evaluate the foregoing effective moduli $\kappa ^{\ast }$ and $\mu ^{\ast }$. Then, recalling the equivalent interfacial relations (2.11) and (2.12), Eqs (3.3) and (3.4) are rewritten based on the divergence theorem as

Accordingly, the strain energy $W$ stored in the RVE model can be computed by

where $\left(\boldsymbol{\bar{\sigma}}\right) ^{\prime }$ and $\left(\boldsymbol{\bar{\varepsilon}}\right) ^{\prime },$ respectively, denote the deviator part of $(\boldsymbol{\bar{\sigma}})$ and $\left(\boldsymbol{\bar{ \varepsilon}}\right)$, and can be calculated by Eq (3.2).

Introducing Eqs (3.1) and (3.6) into Eq (3.7) delivers an intuitive expression of $W$ such that

It is clear from Eq (3.8) that different forms of $\boldsymbol{\varepsilon }^{0}$ can be adopted to determine the overall moduli $\kappa ^{\ast }$ and $\mu ^{\ast }$ of PBXs. Here, $\boldsymbol{\varepsilon }^{0} = \boldsymbol{I}$ and $\left[ 0\text{ }1\text{ } 0;1\text{ }0\text{ }0;0\text{ }0\text{ }0\right]$ are employed and individually prescribed on $\partial \Omega _{RVE}$ to evaluate the averaged values of $\boldsymbol{\bar{\sigma}}$ and $\boldsymbol{\bar{\varepsilon}}$. Hereafter, the effective moduli of PBXs are calculated by

Up to now, the scheme using the micro-mechanics to compute the effective moduli of PBXs has been introduced completely, and will be applied in the following section to deduce the analytical solution of a simplified PBX model containing the imperfect interface model characterized by Eqs (2.11) and (2.12).

4.   A simplified PBX model together with its effective moduli

In this section, we focus on deducing the exact effective moduli of a designated PBX model containing the aforementioned imperfect interface model, which may be served as a benchmark for computational methodology verification and preliminary evaluation of PBX material property. Let's first consider an infinite domain $\Omega$ consisting of a spherical particle $\Omega ^{\left(1\right) }$ of radius $R$ embedded in an infinite matrix domain $\Omega ^{\left(2\right) }$; see Figure 4. Further, $\Omega ^{\left(1\right) }$ and $\Omega ^{\left(2\right) }$ are supposed to be made of linearly elastic explosive material of the same kind, and $\Gamma _{0}$ between $\Omega ^{\left(1\right) }$ and $\Omega ^{\left(2\right) }$ is assumed to be a spring-type interface characterized by the aforementioned relations (2.11) and (2.12), and constituted by the binder material. In addition, the global Cartesian coordinate system $(\boldsymbol{e}_{1}, \boldsymbol{e}_{2}, \boldsymbol{e}_{3})$ is chosen to locate at the center of the spherical particle. Thus, $\Gamma _{0}$ can be described mathematically as

and its unit normal vector $\boldsymbol{n}$ oriented from $\Omega^{(1)}$ into $\Omega^{(2)}$ satisfies

The overall elasticity of the PBX model in Figure 4 is isotropic in view of the symmetry of the geometrical model together with the materials applied. Accordingly, the boundary value problem (BVP) to determine its effective moduli is formulated as follows:

Find the physical fields $\boldsymbol{u}(\boldsymbol{x}), \boldsymbol{ \varepsilon }(\boldsymbol{x})$, and $\boldsymbol{\sigma }(\boldsymbol{x})$ within $\Omega$ fulfilling the conditions such that

($i$) Inside each domain $\Omega ^{\left(r\right) }$ $\left(r = 1, 2\right) :$

($ii$) Relations Eqs (2.11) and (2.12) hold across $\Gamma _{0}$.

($iii$) Two distinct types of boundary conditions, i.e.,

with

are prescribed individually on the remote surface $\partial \Omega$ of the simplified model to evaluate the effective bulk and shear moduli. Above, $\varepsilon _{0}$ is a nonzero constant. In what follows, we shall first derive the analytical responses, i.e., the displacement, strain, and stress distributions, of the foregoing BVPs. Hereafter, the exact solution of the effective moduli of the model in Figure 4 are derived by using the generalized self-consistent scheme.

4.1.   Exact elastic responses of the BVP involving the specified boundary conditions

In what follows, we shall deduce the analytical solutions to the foregoing BVPs governed by Eqs (4.3)–(4.6).

4.1.1.   Analytical solution of the BVP involving the boundary $BC_{1}$

Let's first consider the BVP under the remote uniform strain condition (4.4) together with (4.5). Inspired by ^[33,35], the displacement field in each domain $\Omega ^{\left(r\right) }$ of the benchmark model is postulated as

with $\alpha ^{(r)}$ and $\beta ^{(r)}$ being the four unknowns to be determined. Substituting Eq (4.7) into the aforementioned governing equation (4.3), it is readily shown that the very local equilibrium equations, i.e., ${div}\boldsymbol{\sigma }^{(r)} = \mathbf{0, }$ are satisfied automatically. Hence, these displacement fields are the very solution to the BVP under investigation.

Then, we shall determine the unknown coefficients by using the boundary conditions. Introducing Eq (4.7) of $\Omega ^{\left(2\right) }$ into the remote boundary condition (4.4) leads to $\alpha ^{(2)} = 1$. Further, the displacement field $\boldsymbol{u}^{(1)}$ is physically meaningful only when its value at the center of $\Omega ^{\left(1\right) }$ is finite. We thus obtain $\beta ^{(1)} = 0$.

Equation (4.7) is then greatly simplified by taking account of $\alpha ^{(2)} = 1$ and $\beta ^{(1)} = 0$ as

Next, we need to determine the rest unknowns $\alpha ^{\left(1\right) }$ and $\beta ^{(2)}$ by using the aforementioned interfacial relations (2.11) and (2.12) across $\Gamma _{0}$. On use of Eq (4.3), we can derive the secondary physical quantities, i.e., the strain and stress fields, over $\Omega ^{\left(r\right) }$ as

Combining Eqs (2.11), (2.12), (4.9), and (4.10) yields

Now, the two unknowns can be uniquely determined by solving these two linearly independent equations (4.11) and explicitly expressed as

Thus far, the exact displacement, strain, and stress fields of each domain $\Omega ^{\left(r\right) }$ in Figure 4 are obtained by substituting the coefficients in Eq (4.12) back into Eqs (4.8)–(4.10).

4.1.2.   Analytical solution of the BVP involving the boundary $BC_{2}$

We shall solve the BVP under the remote uniform shear strain boundary condition equation (4.4) combined with (4.6) in this subsection. The displacement field in each phase $\Omega ^{\left(r\right) }$ satisfies

where $a_{1}, b_{1}, b_{2}$, and $c_{2}$ are the model coefficients to be determined. It is easy to check that the displacement fields adopted satisfy the very local equilibrium equations, i.e., ${div}\boldsymbol{\sigma }^{(r)} = \mathbf{0}$. Interested readers are referred to ^[31,32,37] for more details of the deductions. Thus, the foregoing displacement fields are just the solution to the BVP in question.

Likewise, the former subsection, the relevant secondary physical fields, i.e., the strain and stress fields, are also deduced at first and the detailed expressions are presented in the Appendix. Introducing the physical fields (4.13) and (4.14) together with the secondary physical fields into the relations (2.11) and (2.12) across $\Gamma _{0}$ delivers 4 linearly independent equations such that

with $c_{11}$, $c_{12}$, $...$, and $c_{45}$ being the coefficients defined in the Appendix.

Solving the foregoing four equations, we eventually determine the four coefficients as

where $\lambda _{0}$, $\lambda _{11}$, $\lambda _{12}$, $...$, and $\lambda _{34}$ are the coefficients detailed in the Appendix.

Up to now, the exact displacement, strain, and stress fields of each domain $\Omega^{\left(r\right) }$ in Figure 4 under the remote shear strain boundary condition are obtained by substituting the coefficients in Eq (4.16) back into Eqs (4.13) and (4.14).

4.2.   Analytical effective isotropic moduli of the benchmark model

With the aid of the analytical strain/stress solutions under the aforementioned remote boundary conditions, a generic micro-mechanical scheme ^[37,39,40], such as the Mori-Tanaka method, generalized self-consistent method, and so forth, can be applied to evaluate the effective moduli of the foregoing benchmark model. In the current paper, the equivalent scheme established by Duan et al. ^[37] is employed to calculate these quantities. More precisely, two steps necessitate being performed to achieve this goal: (i) replacing the embedded particle coated by the binder interface to be an equivalent inhomogeneity covered by a perfect interface; (ii) computing the effective moduli via the generalized self-consistent scheme established in ^[41]. We shall derive the formula of the benchmark model in Figure 4 containing the spring-type interface equations (2.11) and (2.12) in the following subsections.

4.2.1.   Equivalent replacement based on the energy principle

The primary idea of the equivalent replacement by Duan et al. ^[37] lies in that the elastic energy stored in a particle/imperfect interface/coating system must be equal to that in an equivalent-particle/perfect interface/coating system regardless of the remote uniform boundary condition applied. As for an infinite model consisting of pure coating material under the remote uniform strain condition equation (4.4), the physical fields within $\Omega$ can be trivially written as

Accordingly, the elastic energy $E_{h}$ stored in $\Omega$ can be expressed as

where $V$ refers to the volume of the infinite model $\Omega$.

Next, we cut out a sphere of radius $R$ from $\Omega$; meanwhile, substituting back an equivalent particle $\Omega ^{(1)}$ of radius $R$. Precisely, $\Omega ^{(1)}$ is perfectly bonded to the coating domain $\Omega ^{(2)}$ through the interface $\Gamma _{0}$. The material forming $\Omega ^{(1)}$ is linearly isotropic, thus,

Compared with the foregoing homogeneous model, the elastic energy stored in this particle/perfect interface/coating system equals

with

Above, $V_{1}$ denotes the volume of $\Omega^{(1)}$, and $\mathbb{S}$ is the classical interior Eshelby tensor (see, e.g., ^[40]) defined by

Then, we recall the practical case where $\Omega ^{(1)}$ is bonded to $\Omega ^{(2)}$ through an energetic interface $\Gamma _{0}$ characterized by the relations (2.11) and (2.12). Under the remote uniform strain boundary condition, the elastic energy stored within such a particle/imperfect interface/coating system can be evaluated through ^[35]

with

In the above formula, $\boldsymbol{\sigma }^{\left(2\right) }$ and $\boldsymbol{ u}^{\left(2\right) }$ separately denote the stress and displacement values evaluated at the interface $\Gamma _{0}$, but from the coating $\Omega ^{(2)}$ side.

Combining the energy equivalency condition below,

we can eventually determine the moduli $\kappa_{eq}$ and $\mu _{eq}$ of the material forming $\Omega^{(1)}$ within the equivalent particle/perfect interface/coating system.

4.2.2.   Equivalent isotropic moduli of the benchmark model

To calculate the equivalent bulk modulus $\kappa_{eq}$, we invoke the aforementioned practical BVP under the remote uniform strain condition (4.4) combined with (4.5). The physical fields within $\Omega^{(r)}$ are given by Eqs (4.8)–(4.12). Then, introducing Eqs (4.8), (4.10), and (4.17) into Eq (4.23) delivers

Inserting Eqs (4.5), (4.18), and (4.21) into Eq (4.20) provides

Accounting for the energy equivalency relation equation (4.24) immediately gives the exact expression of $\kappa_{eq}$ as follows:

When $\mu_{eq}$ of the equivalent particle/perfect interface/coating system is concerned, we need to revisit the heterogeneous model in Figure 4 under the remote uniform strain boundary (4.4) with (4.6). In a similar fashion, combining Eqs (4.14), (6.5), (4.17), and the integral relation (4.23), we obtain

Substituting Eqs (4.6), (4.18), and (4.21) into (4.20) brings

Finally, the aforementioned energy equivalency condition Eq (4.24) is applied and results in

Thus far, the two equivalent isotropic moduli of the equivalent particle/perfect interface/coating system have been derived and will be applied to evaluate the effective moduli of the foregoing benchmark model.

4.2.3.   The overall isotropic moduli of the benchmark model

It is now the time to determine the effective moduli $\kappa ^{\ast }$ and $\mu ^{\ast }$ of the designated benchmark model. Among various micormechanics methods ^[39,40], we adopt the generalized self-consistent method (GSCM) proposed in ^[41] to compute these parameters. The effective bulk modulus $\kappa ^{\ast }$ of a heterogeneous model identical to Figure 4 but containing an equivalent particle covered by perfect interface $\Gamma_{0}$ can be calculated by ^[40,41]

where $f$ represents the volume fraction of the spherical particle domain $\Omega ^{\left(1\right) },$ and $\kappa_{eq}$ is given in Eq (4.27).

Referring to the effective shear modulus $\mu ^{\ast }$, its estimation by the generalized self-consistent model (GSCM) is the positive root of the following quadratic equation ^[41]:

with $\eta _{1}$, $\eta _{2}$, and $\eta _{3}$ the model coefficients defined by

Above, ${\nu }_{p}$ is the Poisson's ratio of the material forming $\Omega^{(2)}$, and $\mu ^{\left(eq\right) }$ is expressed in Eq (4.30).

Up to now, the analytical effective isotropic moduli $\kappa ^{\ast }$ and $\mu ^{\ast }$ of the designated PBX model have been explicitly presented in Eqs (4.27), (4.30)–(4.33) together with (4.12)–(4.16).

4.3.   Upper/lower bounds of the effective moduli of the benchmark model

As implied by the PBX model geometry and boundary condition, the foregoing analytical effective isotropic moduli are derived based on the assumption that the model is in infinite dimension and the embedded particle domain $\Omega^{(1)}$ is much smaller than the outer coating domain $\Omega^{(2)}$. Here, the upper and lower bounds of the effective elastic moduli of the designated PBX model are also provided for the purpose of the analytical result validation. Invoking the method initially established by Hashin ^[31] and extended by Zhu et al. ^[36] to determine the upper/lower bounds of effective moduli of matrix-based composites, we can derive the explicit expression of the bounds, namely, $\overline{\kappa }^{\ast }$, $\underline{\kappa }^{\ast }$ and $\overline{\mu }^{\ast }, \underline{\mu }^{\ast },$ of the aforementioned PBX model containing spring-type interface $\Gamma_{0}$ as

Above, $f$ also refers to the volume fraction of the particle domain $\Omega ^{\left(1\right) }$. Interfacial parameters $\omega _{n}$ and $\omega _{T}$ are defined in Eq (2.13).

5.   Verification and discussion

This section is dedicated to examining the correctness of the derived exact effective moduli of the foregoing simplified PBX model, and discussing the influence of particulate size on the overall elasticity of PBXs.

5.1.   Analytical solution examination

Let's revisit the PBX model in Figure 4, and consider it as an RVE of PBXs of a certain kind. The side length of $\Omega$ is set as $w = 8\times 10^{-4}m$ and the radius $R$ of the spherical particle $\Omega^{(1)}$ takes the value ranging from $0.01 w$ to $0.4\; w$, since the characteristic dimension of practical explosive particles varies from several to hundreds of micrometers (c.f., ^{[1,3,4,6,9,42]}). Besides, two different situations separately with $f_{p} = 90\%$ and $95\%$ are taken as instances for the afterward verification. To be noted, $f_{p}$ refers to the volume fraction of the explosive particles, and approximately equals

Above, $h$ denotes the thickness of the equivalent binder interface $\Gamma_{0}$.

Domains $\Omega^{(1)}$ and $\Omega^{(2)}$ are supposed to be constituted by the commonly applied octogen(HMX) crystal of the same kind, and its elastic parameters (see ^[6]) are listed in Table 1. To cover a large portion of binder materials of practical interest ^[6,42], the elastic moduli of the substance forming $\Gamma_{0}$ are selected as $E_{b} = (10^{-7}\sim10^{-1})E_{p}$ and $v_{b} = 0.15$.

Figures 5 and 6 plot the analytical evaluations and numerical predictions of the effective isotropic moduli together with their upper and lower bounds of the aforementioned PBX model. It is obvious that the analytic exact predictions of each instance by Eqs (4.7)–(4.33) are rigorously located between the relevant upper and lower bounds. As the radius $R$ of $\Omega^{(1)}$ decreases, as expected, the effective bulk (or shear) modulus $\kappa^{*}$ (or $\mu^{*}$) gradually approaches the bulk (or shear) modulus $\kappa_{p}$ (or $\mu_{p}$) of the explosive particles. Furthermore, the numerical predictions by the XFEM fit pretty well with the analytic approximations for all combinations of the binder and explosive particle material as the radius $R$ of the embedded particle is smaller than $0.3\; w$. When $R$ get larger, the analytical solutions gradually deviate from the numerical predictions due to fact that they are derived based on the hypothesis that the inclusion volume fraction of the embedded explosive particle is much smaller than the volume fraction of the outer explosive domain. It is thus concluded that the derived analytic exact expressions are correct, and the foregoing equivalent replacement strategy combined with the GSCM method is efficient to predict the effective moduli of PBXs.

To be noted, the overall moduli of PBXs are significantly $E_{b}$ dependent, particularly when the elastic moduli of the binder and explosive particles are in high contrast, say $E_{b} = 10^{-7}E_{p}$. More significant degradation of the effective moduli $\kappa^{*}$ and $\mu^{*}$ shows up as $E_{b}$ decreases according to Figures 5 and 6. Hence, it's of particular importance to take into account the quantitative influence of the binder properties as PBXs are designed for practical applications. It is suggested that the binder material with smaller $E_{b}$ values be selected so that the binder phase may generate large deflections during the fabrication and transportation process, which protect the explosive particles from serious damages.

In addition, the upper/lower bounds of $\kappa^{*}$ (or $\mu^{*}$) cover a very wide range as $E_{b}$ is lower than $E_{p}$ by 4 or more orders of magnitude; see Figures 5(a), (b) and 6(a), (b). It is pretty hard to roughly evaluate the PBX moduli information by using the bound formula under such situations. When $E_{b}$ is close to $E_{p}$, (see Figures 5(c), (d) and 6(c), (d)), variations of the PBX moduli between the upper and lower bounds get pretty tight, particularly as the radius $R$ of the embedded particle is larger than $0.4\; w$. Thus, the two bounds under this situation are suitable to predict the effective moduli of PBXs in practical applications.

5.2.   Size-dependent effective moduli of PBXs

We now extend the simplified PBX model in Figure 4 to further investigate the size-dependent effect of PBXs. The side length $w$ of the model this time is assumed to vary from $5\times10^{-8}$ m to $8\times10^{-4}$ m, and the radius $R$ of the inner particle $\Omega^{(1)}$ fulfills $R = 0.35\; w$. The thickness $h$ of the binder interface $\Gamma_{0}$ is considered as being uniform and equals $10^{-8}$ m.

Both the explosive particle domains $\Omega^{(1)}$ and $\Omega^{(2)}$ are supposed to be made of the foregoing HMX crystals, and the relevant elastic coefficients are listed in Table 1. Likewise to the former subsection, the elastic moduli of the binder material forming $\Gamma_{0}$ satisfy the relations such that $E_{b} = (10^{-7}\sim0.1)E_{p}$ and $v_{b} = 0.15$.

The effective moduli $\kappa^{*}$ and $\mu^{*}$ by Eqs (4.7)–(4.33) against $w$ are plotted in Figure 7. It is observed that these overall moduli are closely related not only to the size of the explosive particles but also to the stiffness of the binder materials. As $E_{b}$ is smaller than $E_{p}$ by 7 or more orders of magnitude, $\kappa^{*}$ and $\mu^{*}$ become insensitive to the size of the explosive particles for the whole range investigated. With the increment of $E_{b}$, obvious influence induced by the particle size is clearly observed. Precisely, larger $E_{b}$ value makes $\kappa^{*}$ and $\mu^{*}$ be more sensitive to smaller explosive particles. All in all, both the binder properties and explosive particle size should be considered when the overall moduli of PBXs are investigated.

5.3.   Influence of the volume fraction of the binder medium on the effective moduli of PBXs

In this subsection, the simplified PBX model in Figure 4 is also applied to examine the influence of the binder volume fraction on the overall moduli of PBXs. The side length $w$ of the model is specified as $10^{-5}$ m, and the radius $R$ of the inner particle $\Omega^{(1)}$ is set as $R = 0.3\; w$. The volume fraction of the binder phase changes in the range $[0.01, 0.1]$, and the thickness $h$ of the binder interface $\Gamma_{0}$ is supposed to be uniform and evaluated by the formula (5.1).

Similar to the former subsection, both the explosive particle domains $\Omega^{(1)}$ and $\Omega^{(2)}$ are considered as being made of the HMX crystals, and the relevant elastic coefficients are listed in Table 1. The elastic moduli of the binder material forming $\Gamma_{0}$ satisfy the relations such that $E_{b} = (10^{-7}\sim0.1)E_{p}$ and $v_{b} = 0.15$.

The effective moduli $\kappa^{*}$ and $\mu^{*}$ by Eqs (4.7)–(4.33) against $f_{b}$ are visualized in Figure 8. It is seen that both $\kappa^{*}$ and $\mu^{*}$ show quite limited influence by the binder volume fraction $f_{b}$ when $E_{b}$ is smaller than $E_{p}$ by 4 or more orders of magnitude. As $E_{b}$ increases to $10^{-2}E_{p}$ or larger, $\kappa^{*}$ and $\mu^{*}$ become sensitive to $f_{b}$ due to the strengthened bonding effect from the interface. Further, $\kappa^{*}$ and $\mu^{*}$ gradually decrease with the increment of $f_{b}$ due to the fact that $E_{b}$ is smaller than $E_{p}$.

5.4.   Influence of the explosive crystals size on the effective moduli of PBXs

This subsection is dedicated to investigating the explosive crystals size effect on the effective moduli of PBXs using the simplified model in Figure 4 as well. The side length $w$ of the model is specified as $10^{-5}$ m, and the radius $R$ of the inner particle $\Omega^{(1)}$ fulfills $R = [0.1, 0.2, 0.3, 0.4]\; w$. The volume fraction of the binder phase remains unchanged $0.05$, and thickness $h$ of the binder interface $\Gamma_{0}$ is supposed to be uniform and evaluated by the formula (5.1).

Similar to the former subsection, both the explosive particle domains $\Omega^{(1)}$ and $\Omega^{(2)}$ are supposed to be made of the foregoing HMX crystals, and the relevant elastic coefficients are listed in Table 1. The elastic moduli of the binder material forming $\Gamma_{0}$ satisfy the relations such that $E_{b} = (10^{-7}\sim0.1)E_{p}$ and $v_{b} = 0.15$.

The effective moduli $\kappa^{*}$ and $\mu^{*}$ by Eqs (4.7)–(4.33) against $f_{b}$ are visualized in Figure 9. We can observe that the overall moduli $\kappa^{*}$ and $\mu^{*}$ are closely related to the size of the embedded explosive particle $\Omega^{(1)}$. With the increment of the radius $R$, both $\kappa^{*}$ and $\mu^{*}$ decreases dramatically, particularly, when the binder Young's modulus $E_{b}$ takes the smaller values.

6.   Concluding remarks

In the present paper, a physically sound spring-type interface model is derived for the first time to characterize the thin binding interphases inside PBXs based on some idealized hypothesis, and the averaged and effective quantities over a RVE of PBXs are reformulated according to the primitive definitions in micro-mechanics. Hereafter, a simplified PBX model is designated, and its exact effective isotropic moduli are obtained on use of the equivalent replacement strategy ^[37] and generalized self-consistent scheme ^[41]. The upper and lower bounds of the designed PBX model are also deduced and given explicitly. Taking the HMX-based PBX as an illustration, we eventually verified the correctness of the exact analytic solution of the effective isotropic moduli. To be specified, the explicit expressions of the analytic effective isotropic moduli and their upper/lower bounds provide a direct tool to estimate the basic elastic properties of PBXs, and can be applied as benchmarks to test the validity of numerical schemes to predict the effective moduli of highly filled composites as well.

Besides, the size-dependent effect of PBXs is also studied with the aid of a simplified PBX model. It is found that the effective isotropic moduli are strongly dependent not only on the size of explosive particles but also on the stiffness of binder material. Therefore, it is a better choice to synthesize practical PBXs by taking into account these two factors simultaneously. Remark that our attention of this study is mainly on the overall moduli evaluation of the PBXs consisting of the energetic crystals showing the isotropic properties. The PBXs involving the anisotropic energetic crystals together with their influence factors will be investigated in our future works.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work is funded by the National Natural Science Foundation of China (Grant numbers 12372142 and 51305362), and the Key Research and Development Project in Sichuan Province (Grant number 2025YFHZ0040).

Conflict of interest

The authors declare there are no conflicts of interest.

Appendix

1) Expression of the secondary physical fields of the BVP in Section 4.1.2.

2) Expression of the coefficients in Eq (4.15).

3) Expression of the coefficients in Eq (4.16).