2.3.1. The Method Solving Strain (εr, εθ)

The differential relation of strain (εr, εθ) and the particle velocity is established by the conservation equation. Now, the initial condition t = 0, v(ri,t)=0 is known, and v(ri,t) at different positions ri (i = 1,2, ... ,) is also known. So, the time numerical integration operation can be performed to determine εθ(ri, t). Then, the first derivative *∂*v(ri,t)/*∂*t can be obtained by numerical differential operation. Similarly, the strain εr(ri,t) can be determined by integrating time.

## 2.3.2. The Method Solving Strain (σr, σθ)

However, the stresses σ<sup>r</sup> and σ<sup>θ</sup> are still unknown. The system composed of volume and shape deformation is to be determined. The solving of σ<sup>r</sup> and σ<sup>θ</sup> in this way is not sufficient, and one of the equation relations must be known first. In the elastic range, it is accepted that the volume deformation satisfies the linear law of elasticity (2a) and is independent of the rate. Then, it is easy to determine this relationship under quasi-static conditions. The calculation process related to quantity ε<sup>r</sup> and εθ, *∂*σθ/*∂*t, and *∂*σr/*∂*t can be expressed in Equations (4) and (5b).

In order to establish the magnitude relationship at each radius, the path-line processing method can be used to define the total derivative of a certain magnitude on the path-line (Grady, 1973), and the path-lines P1, P2, P3 ... Pi ... Pm can be established as shown in the figure. In this way, when the spherical particle velocity histories v(ri,t) at multiple Lagrangian radii r=ri are provided, and their related other time and position differential components *∂*v(ri,t)/*∂*t can be easily determined.

$$\frac{\partial \sigma\_{\theta}}{\partial t} = \frac{1}{2} \left( \left( \Im K \frac{\partial \varepsilon\_{r}}{\partial t} + 2 \frac{\partial \varepsilon\_{\theta}}{\partial t} \right) - \frac{\partial \sigma\_{r}}{\partial t} \right) \tag{4}$$

$$\left.\frac{d\sigma\_r}{dr}\right|\_p = \left.\frac{\partial\sigma\_r}{\partial r}\right|\_t + \left.\frac{\partial\sigma\_r}{\partial t}\right|\_r \frac{dt}{dr} = \left.\frac{\partial\sigma\_r}{\partial r}\right|\_t + \left.\frac{\partial\sigma\_r}{\partial t}\right|\_r \frac{1}{r'}\Big|\_p \tag{5a}$$

substituting (4) into (5a), the calculation formula of partial derivative about stress *∂*σr/*∂*t can be expressed as (5b).

$$\frac{\partial \sigma\_r}{\partial t} = r' \left( \left. \frac{d\sigma\_r}{dr} \right|\_p - \rho\_0 \frac{\partial v}{\partial t} + \frac{2(\sigma\_r - \sigma\_\theta)}{r} \right) \tag{5b}$$

The zero initial condition is known at different positions of wave propagation (σ<sup>r</sup> = 0 along path-line P1), and the stress σ<sup>r</sup> at different radius r=ri along the path-line P2 (Figure 3) is obtained through the integration of partial derivative *∂*σr(ri,t)/*∂*t by using the constructed path-lines (5). Then, *∂*σθ(ri,t)/*∂*t is known from (4), and the circumferential stress at different positions r = ri on the path-line P1 σθ(ri,t)|P=j can be calculated by integrating *∂*σθ(ri,t)/*∂*t. Similarly, the stress σr(ri,t)|P=j+1 and σθ(ri,t)|P=j+1 on all path-lines can be determined by cycling in sequence. Note that this method can be used to load the whole process, which is called "NV + T0/SW" for short.

**Figure 3.** Results of velocity histories in concrete: (**a**) schematic diagram of test location layout in Mid-plane; (**b**) the series of measured particle velocity profiles.
