2.3.2. Algebraic Identifier with Variable Operation Velocity

In this section, the rotor-bearing system velocity is considered as a linear ramp of excitation. The mathematical model of the system is defined by Equation (3). In order to develop the parameter identifier, this equation is rewritten as follows

$$\left[\mathbf{M}\right]\left\{\ddot{\delta}\right\} + \left[\mathbf{C}\_1 + \dot{\phi}\mathbf{C}\_2\right]\left\{\dot{\delta}\right\} + \left[K\_1 + K\_2 + \ddot{\phi}K\_3\right]\{\delta\} = \dot{\phi}^2 F\_1(\phi) + \ddot{\phi}F\_2(\phi) \tag{16}$$

By multiplying Equation (16) by *t* 2 and integrating the result twice with respect to time, the following is obtained

$$\int\_{}^{\langle 2 \rangle} \left[ [M] \left\{ \bar{\delta} \right\} + \left[ \mathbb{C}\_1 + \dot{\phi} \mathbb{C}\_2 \right] \left\{ \dot{\delta} \right\} + \left[ K\_1 + K\_2 + \bar{\phi} K\_3 \right] \{ \delta \} \right] t^2 = \int\_{}^{\langle 2 \rangle} \left\{ \dot{\phi}^2 \mathbb{F}\_1(\phi) + \bar{\phi} \mathbb{F}\_2(\phi) \right\} t^2 \tag{17}$$

where 0 (2) *ϕ*(*t*) are iterated time-integrals of the form 0 *t* 0 0 *σ*1 0 ··· 0 *σn*−1 0 *ϕ*(*σn*)*dσn* ··· *dσ*1 with 0 *ϕ*(*t*) = 0 *t* 0 *ϕ*(*σ*)*dσ*, and *n* a positive integer.

Similarly for the case of constant velocity, matrices [*K*1] and [*C*1] contain the stiffness and damping parameters provided by the supports. Therefore, after the integration of the left part of Equation (17) and rearranging terms, we have

$$\begin{aligned} &\int \prescript{(2)}{}{\{\mathcal{K}\}}\_{1} t^{2} \{\delta\} + [\mathcal{C}\_{1}] \left[ \int t^{2} \{\delta\} - 2 \int ^{(2)}{}{t} \{\delta\} \right] = \int \left[ 4Mt - \dot{\phi} \mathcal{C}\_{2} t^{2} \right] \{\delta\} \\ &+ \int \prescript{(2)}{}{\{\mathcal{C}\_{2} \left(\ddot{\phi} t^{2} + 2\dot{\phi} t\right) - 2M - \left(K\_{2} + \ddot{\phi} K\_{3}\right) t^{2}\}} \{\delta\} - [M] t^{2} \{\delta\} \\ &+ \int \prescript{(2)}{}{\{\dot{\phi}^{2} F\_{1}(\phi) + \ddot{\phi} F\_{2}(\phi)\}} t^{2} \end{aligned} \tag{18}$$

It is worth mentioning that Equation (18) can be separated into individual equation systems for each node where the bearings are located. These equations can be written as follows

$$
\begin{bmatrix} k\_{xx} & k\_{xz} \\ k\_{zx} & k\_{zx} \end{bmatrix} \int^{\langle 2 \rangle} t^2 \begin{Bmatrix} u\_i \\ & w\_i \end{Bmatrix} + \begin{Bmatrix} c\_{xx} & c\_{xz} \\ c\_{zx} & c\_{xz} \end{Bmatrix} \left( \int t^2 \begin{Bmatrix} u\_i \\ & w\_i \end{Bmatrix} - 2 \int^{\langle 2 \rangle} t \begin{Bmatrix} u\_i \\ & w\_i \end{Bmatrix} \right) = \begin{Bmatrix} b\_{ni} \\ b\_{ni} \end{Bmatrix} \tag{19}
$$

To solve Equation (19), an equal number of equations and unknows is required. Therefore, Equation (19) is successively integrated three times to obtain the missing equations which are expressed as

$$
\begin{aligned}
\begin{bmatrix} k\_{xx} & k\_{xz} \\ k\_{zx} & k\_{zz} \end{bmatrix} \end{aligned}
\begin{aligned}
\begin{aligned}
\begin{Bmatrix} (3) \\ w\_{i} \end{Bmatrix} \end{aligned} \\
\begin{aligned}
\begin{Bmatrix} \boldsymbol{w}\_{i} \\ \boldsymbol{w}\_{i} \end{Bmatrix} \end{aligned} + \begin{bmatrix} \mathbf{c}\_{xx} & \mathbf{c}\_{xz} \\ \mathbf{c}\_{xx} & \mathbf{c}\_{zz} \end{bmatrix} \left( \int \begin{array}{c} \mathbf{t}^{(2)} \\ \mathbf{w}\_{i} \end{array} \begin{Bmatrix} \mathbf{u}\_{i} \\ \mathbf{w}\_{i} \end{Bmatrix} \right) - 2 \int \begin{Bmatrix} \mathbf{u}\_{i} \\ \mathbf{w}\_{i} \end{Bmatrix} \end{aligned} \\
\begin{aligned}
\begin{Bmatrix} \mathbf{u}\_{i} \\ \mathbf{w}\_{i} \end{Bmatrix} \end{aligned} \tag{20}
\end{aligned}
\end{aligned}
\end{aligned}
\begin{aligned}
\begin{aligned}
\begin{Bmatrix} \mathbf{w}\_{i} \\ \mathbf{w}\_{i} \end{Bmatrix} \end{aligned}$$

$$
\begin{bmatrix} k\_{\text{xx}} & k\_{\text{xz}} \\ k\_{\text{zx}} & k\_{\text{zz}} \end{bmatrix} \int^{(4)} \boldsymbol{t}^2 \left\{ \begin{array}{c} \boldsymbol{w}\_{\text{i}} \\ \boldsymbol{w}\_{\text{i}} \end{array} \right\} + \begin{bmatrix} \boldsymbol{c}\_{\text{xx}} & \boldsymbol{c}\_{\text{xz}} \\ \boldsymbol{c}\_{\text{zx}} & \boldsymbol{c}\_{\text{zz}} \end{bmatrix} \left( \int^{(3)} \boldsymbol{t}^2 \left\{ \begin{array}{c} \boldsymbol{w}\_{\text{i}} \\ \boldsymbol{w}\_{\text{i}} \end{array} \right\} - 2 \int^{(4)} \boldsymbol{t} \left\{ \begin{array}{c} \boldsymbol{w}\_{\text{i}} \\ \boldsymbol{w}\_{\text{i}} \end{array} \right\} \right) = \int^{(2)} \left\{ \begin{array}{c} \boldsymbol{b}\_{\text{ui}} \\ \boldsymbol{b}\_{\text{ui}} \end{array} \right\} \tag{21}
$$

$$
\begin{bmatrix} k\_{xx} & k\_{xz} \\ k\_{zx} & k\_{zz} \end{bmatrix} \int \stackrel{(5)}{\cdot} t^2 \left\{ \begin{array}{cc} u\_i \\ w\_i \end{array} \right\} + \begin{bmatrix} c\_{xx} & c\_{xz} \\ c\_{zx} & c\_{zz} \end{bmatrix} \left( \int \stackrel{(4)}{\cdot} t^2 \left\{ \begin{array}{cc} u\_i \\ w\_i \end{array} \right\} - 2 \int ^{(5)} t \left\{ \begin{array}{cc} u\_i \\ w\_i \end{array} \right\} \right) = \int ^{(3)} \left\{ \begin{array}{cc} b\_{\text{ini}} \\ b\_{\text{ini}} \end{array} \right\} \tag{22}
$$

From Equations (19)–(22), a linear system equation is obtained for each node where the bearings are located. These equations can be expressed as

$$[A\_{\mathfrak{s}}(t)]\{\Theta\_{\mathfrak{s}}\} = \{b\_{\mathfrak{s}}(t)\}\tag{23}$$

where {Θ*s*} = {*kxx kxz kzx kzz cxx cxz czx czz* }*T* denotes the transposed vector of parameters to be identified and [*As*(*t*)], {*bs*(*t*)} are 8 × 8 and 8 × 1, respectively.

Again, the condition *det*[*As*(*t*)] = 0 must be satisfied to identify the vector {Θ*s*}.

From the solution of Equation (23), a mathematical model for an on-line identifier of stiffness and damping bearing parameters can be obtained as

$$\{\Theta\_{\mathfrak{s}}\} = \left[A\_{\mathfrak{s}}^{-1}\right] \{b\_{\mathfrak{s}}\} \,\,\forall t \in (t\_{0}, t\_{0} + \mathfrak{e}] \tag{24}$$

As can be observed, algebraic identification of stiffness and damping bearing parameters is independent of system initial conditions and only depends on the displacement vector and the type of ramp excitation. It is important to mention that as with the case of constant velocity, to identify the supports parameters, lateral vibration measurements at the node and the nodal slopes are required. Moreover, similar information from the adjacent node is also needed.

## **3. Results**

In Figure 3, a scheme of the rotor-bearing system considered in this work and its discretization is presented.

**Figure 3.** Rotor-bearing system scheme [40].

To obtain the mathematical model for the rotor-bearing system using the FE method, it was discretized into 11 beam-like elements, as is shown in Figure 3. The system includes two inertial disks located at nodes 1 and 12, while supports (bearings) are placed at nodes 4 and 8. The correct nodal location ensures that the simulation replicates the model's real conditions. In addition, two unbalanced masses were considered in two different angular positions located on inertial disks *D*1 and *D*2.

In Table 1, the mechanical and geometrical properties of the shaft are shown, while the inertial properties of discs and unbalanced masses are presented in Table 2.

**Table 1.** Mechanical and geometrical properties of the rotor-bearing shaft.

**Table 2.** Inertial properties of the disks and unbalance masses.

In Table 3, the stiffness and damping bearing parameters [40] used for numerical simulation are presented.

**Table 3.** Stiffness and damping bearing parameters [40].

On-line algebraic identification of stiffness and damping bearing parameters was determined based on the vibratory response of the rotor-bearing system in the time domain, which was obtained from Equations (3) and (7) by using the Newmark method for numerical integration.

#### *3.1. Algebraic Parameter Identification with Constant System Velocity*

The displacement vector used in the algebraic identification procedure was obtained from Equation (7) by using the Newmark method for numerical integration and taking into account a constant rotational velocity of the rotor-bearing system.

In Figure 4, vibration signals at node 4 (corresponding to bearing 1 location) of the rotor-bearing system of Figure 3 are presented. This response is obtained for an operation rotational velocity Ω = 600 rpm. These signals, the nodal slopes and the corresponding information of the nodes adjacent to the bearing locations are the required data to identify stiffness and damping parameters of the bearings.

**Figure 4.** Vibration signal at node 4 (bearing 1) at 600 rpm: (**a**) X direction; (**b**) Z direction.

Figures 5 and 6 present the obtained results from the numerical simulation for the algebraic identification of stiffness and damping parameters for bearing 1, while Figures 7 and 8

show the results corresponding to bearing 2. It is important to mention that the sample time used in the simulation was 0.1 milliseconds. However, by carrying out numerical simulations with different sample times, it was observed that the shorter the sampling period, the faster the identifier converges.

**Figure 5.** Identified stiffness parameters for bearing 1 at 600 rpm. (**a**) *kxx*, (**b**) *kxz*, (**c**) *kzx*, (**d**) *kzz*.

**Figure 6.** Identified damping parameters for bearing 1 at 600 rpm. (**a**) *cxx*, (**b**) *cxz*, (**c**) *czx*, (**d**) *czz*.

**Figure 7.** Identified stiffness parameters for bearing 2 at 600 rpm. (**a**) *kxx*, (**b**) *kxz*, (**c**) *kzx*, (**d**) *kzz*.

**Figure 8.** Identified damping parameters for bearing 2 at 600 rpm. (**a**) *cxx*, (**b**) *cxz*, (**c**) *czx*, (**d**) *czz*.

As can be observed in Figures 5–8, the identification of both stiffness and damping parameters of the bearings is carried out in less than 0.1 s, and once the parameter reaches the identified value, this remains for the rest of the time period. For a better analysis of the identifier behavior, only results for 0.1 s are presented in Figures 5–8, because it is important to observe the time that the identifier requires to converge to the estimated value.

### *3.2. Algebraic Parameter Identification with Variable System Velocity*

The displacement vector used as input data for the algebraic identification is obtained from Equation (3) by using the Newmark method for numerical integration and taking into account a linear ramp excitation with angular acceleration .. *φ* = 10 *rad*/*s*2. The rotor-bearing system response at node 4 is shown in Figure 9 where the vibratory behavior of the system in the location of bearing 1 can be appreciated.

**Figure 9.** System vibratory response at node 4 with a linear ramp of excitation of 10 *rad*/*s*2.

In Figures 10–13, the behavior of the algebraic identifier for bearing stiffness and damping parameters of both bearings (placed at nodes 4 and 8) is shown as a function of time.

**Figure 10.** Identified stiffness parameters for bearing 1 with a linear ramp of excitation of 10 *rad*/*s*2. (**a**) *kxx*, (**b**) *kxz*, (**c**) *kzx*, (**d**) *kzz*.

**Figure 11.** Identified damping parameters for bearing 1 with a linear ramp of excitation of 10 *rad*/*s*2. (**a**) *cxx*, (**b**) *cxz*, (**c**) *czx*, (**d**) *czz*.

**Figure 12.** Identified stiffness parameters for bearing 2 with a linear ramp of excitation of 10 *rad*/*s*2. (**a**) *kxx*, (**b**) *kxz*, (**c**) *kzx*, (**d**) *kzz*.

**Figure 13.** Identified damping parameters for bearing 2 with a linear ramp of excitation of 10 *rad*/*s*2. (**a**) *cxx*, (**b**) *cxz*, (**c**) *czx*, (**d**) *czz*.

As in the case of constant rotor system velocity, stiffness and damping bearing parameters are identified in les tan 0.1 s as can be observed in Figures 10–13. Furthermore, the parameter values remain constant until the rotor-bearing system reaches its nominal operation velocity. For a better analysis of the identifier behavior, the results for 0.1 s are presented in Figures 10–13 because it is important to observe the required time for the identifier convergence. The sample time used to solve Equation (3) using the Newmark method was 0.1 milliseconds. The numerical solution of Equation (3) was used as input data for the proposed algebraic identifier. Moreover, achieving this sample time with diverse data acquisition systems for experimental implementation was verified.

#### **4. Discussion**

Different numerical simulations were carried out in order to determine the robustness of the proposed identifiers under different conditions for the rotor-bearing system velocity. For the constant velocity case, different magnitudes for the rotor system velocity were considered, while for the variable velocity case, different ramps of excitation were explored.

Figure 14 shows results for the algebraic identification of stiffness and damping parameters for bearing 1 at a constant operation velocity of the rotor-bearing system of 50,000 rpm. A rapid identifier convergence to the estimated values can be observed, meaning that an increase in operation velocity does not affect the identifier performance. It is important to mention that, while the results for the identification of damping parameters of bearing 1 and the stiffness and damping parameters of bearing 2 are not presented, these parameters are correctly identified in less than 0.1 s.

The identifier performance for different ramps' excitation was analyzed. The acceleration values considered for numerical simulation were: .. *φ* = 10 *rad*/*s*2, .. *φ* = 100 *rad*/*s*2, .. *φ* = 1000 *rad*/*s*2, .. *φ* = 3000 *rad*/*s*2 and .. *φ* = 6000 *rad*/*s*2. The result for .. *φ* = 10 *rad*/*s*2 were reported in the previous section. Due to the similar behavior of the identifier with different acceleration values, only results for .. *φ* = 6000 *rad*/*s*2 are shown here. The rotor-bearing system response for a ramp of excitation with the mentioned value of acceleration at node 4 (bearing 1 location) is presented in Figure 15. It can be seen that there is a considerable change in the time scale in comparison with Figure 9 because the acceleration is increased 600 times.

**Figure 14.** Identified stiffness parameters for bearing 1. (**a**) *kxx*, (**b**) *kxz*, (**c**) *kzx*, (**d**) *kzz* at 50,000 rpm.

**Figure 15.** System vibratory response at node 4 with a linear ramp excitation with acceleration of 6000 *rad*/*s*2.

The algebraic identification performance under the conditions described above is shown in Figures 16–19 where the estimation for the stiffness and damping bearings parameter is visualized. For this simulation, the system response in Figure 15 is used as input data.

**Figure 16.** Identified stiffness parameters for bearing 1 with a ramp excitation of 6000 *rad*/*s*2. (**a**) *kxx*, (**b**) *kxz*, (**c**) *kzx*, (d) *kzz*.

**Figure 17.** Identified damping parameters for bearing 1 with a ramp excitation of 6000 *rad*/*s*2. (**a**) *cxx*, (**b**) *cxz*, (**c**) *czx*, (**d**) *czz*.

**Figure 18.** Identified stiffness parameters for bearing 2 with a ramp excitation of 6000 *rad*/*s*2. (**a**) *kxx*, (**b**) *kxz*, (**c**) *kzx*, (d) *kzz*.

**Figure 19.** Identified damping parameters for bearing 2 with a ramp excitation of 6000 *rad*/*s*2. (**a**) *cxx*, (**b**) *cxz*, (**c**) *czx*, (**d**) *czz*.

From the results presented in Figures 16–19, it can be observed that despite the linear ramp excitation being 600 times faster than the corresponding one in Figure 9, the proposed identifier rapidly converges to the estimated values and remains in these values for the rest of the time period. Note that the algebraic identifier is not affected by the system's acceleration and only depends on the displacement vector at each instance. The robustness of the algebraic identification method against acceleration ramp variations had already been proved by Mendoza-Larios et al. [36] but only for the identification of unbalance parameters in rotor-bearing parameters.

Furthermore, the obtained results for both cases, constant and variable rotor-bearing velocity, show a transient state of the identifiers before the convergence to the estimated values of stiffness and damping bearing parameters. This behavior is due to the sample time used in the numerical simulations of the identifiers for solving the iterated integrals of Equations (15) and (24), which utilize the trapezoidal rule. According to Kharab and Guenther [41], this method presents major calculation errors in comparison with other integration methods. However, it was found that the smaller the sample time the shorter the error in the trapezoidal rule calculation.

#### **5. Conclusions**

In this article the identification problem for stiffness and damping parameters of the supports in rotor-bearing systems was addressed. The system model was obtained by the finite element method and using a finite element type beam, which consider the effects of rotational inertia, gyroscopic moments and shear deformations. The algebraic identification technique was applied to the finite element model to obtain two identifiers for the stiffness and damping parameters attributable to the bearings. The first identifier considers a rotor-bearing system operating at a constant velocity, and the second with a linear ramp of excitation as a system velocity input. The numerical results present the identifier behavior showing a fast convergence and robustness in both operation conditions with different values of constant rotational velocity and ramp of acceleration. The numerical results indicate a fast convergence in the stiffness and damping parameters identification in less than 0.06 s for both considered operation conditions. It is important to mention that the convergence time of the identifier depends mainly on the sample time used in numerical simulations. An important characteristic of the proposed algebraic identifiers is that the unbalance parameters (magnitude and phase) are not needed for their development and implementation because only the vibratory response of the system at the bearings' location and the adjacent nodes is required. As a first approach we have proved the proposed identifiers in rotor-bearing system models with constant rotordynamic coefficients. However, as a future work, the proposed identifiers can be used to numerically and experimentally determine rotordynamic coefficients, which are a function of the system rotational velocity as in the case of pressurized bearings, by adapting the identifier method to estimate non-constant functions. This is possible because in the mathematical model used for the identifiers' development, only the effects (stiffness and damping) of the supports are considered without taking into account the nature of the bearings.

**Author Contributions:** Conceptualization, J.G.M.-L., M.A.-M. and S.J.L.-D.; methodology, J.G.M.-L., E.B. and J.C.-O.; software, J.G.M.-L. and S.J.L.-D.; validation, J.G.M.-L., S.J.L.-D. and L.A.B.-T.; formal analysis, J.G.M.-L., M.A.-M. and R.T.-H.; investigation, J.G.M.-L., S.J.L.-D. and L.A.B.-T.; resources, J.G.M.-L., E.B., M.A.-M., R.T.-H. and J.C.-O.; data curation, J.G.M.-L., S.J.L.-D. and L.A.B.-T.; writing original draft preparation, J.G.M.-L., E.B. and M.A.-M.; writing—review and editing, J.G.M.-L., E.B. and M.A.-M.; supervision, E.B. and M.A.-M.; project administration, J.G.M.-L.; funding acquisition, J.G.M.-L. and M.A.-M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** The datasets generated and supporting the findings in the article are obtainable from the corresponding author upon reasonable request.

**Acknowledgments:** The APC was funded by PRODEP-SEP, Mexico.

**Conflicts of Interest:** The authors declare no conflict of interest.


**6. References**

<sup>\*</sup> Corresponding Author

In this chapter, physical phenomena of a MR fluid damper have been carefully investigated through both experimental data and modeling methodologies. Furthermore, a novel forcesensorless control method was proposed and successfully applied to the damping system using the MR fluid damper and without using any mechanical force sensor. The designed force control system is based on the optimized BBM model, functioning as the virtual force

The two test rigs using the MR fluid damper have been fabricated in order to design the models as well as to evaluate the proposed controller. The BBM model is built as the simple direct modeling method and optimized by using neural network technique. The IBBM controller is the direct inverse model of the optimized BBM. In addition, the IBBM has the online self-learning ability by using neural network, consequently, improving the system performance. As the result, the optimized BBM and IBBM models have a strong ability to be applied to damping systems using the MR fluid damper as a virtual force sensor and a

The modeling as well as control validation process has been done by both the simulations and real-time experiments. Based on the simulation and experimental results, it is clear that the optimized BBM can predict the force-displacement behavior of the MR damper with high precision while the IBBM can control the damping system to follow the desired damping force well. As a result, the proposed force-sensorless damping control methodology can become an effective and economical solution for practical MR fluid

[1] Karnopp D, Crosby M.J, Farwood R.A (1974) Vibration control using semi-active force

[2] Yi K, Song B.S (1999) A new adaptive sky-hook control of vehicle semi-active

[3] Kawashima K, Unjoh S, Shimizu K (1992) Experiments on Dynamics Characteristics of Variable Damper. Proc. of the japan national symp. on structural response control,

[4] Mizuno T, Kobori T, Hirai J, Matsunaga Y, Niwa N (1992) Development of Adjustable Hydraulic Dampers for Seismic Response Control of Large Structure. ASME PVP conf.:

suspensions. Proc. of IMechE, Part D: j. auto. eng. 213(3): 293-303.

sensor, and the IBBM model, functioning as the self-learning force controller.

**5. Conclusions**

simple direct controller, respectively.

damper applications in the near future.

*School of Mechanical Engineering, University of Ulsan, Korea*

generators. ASME j. eng. ind. 96(2): 619-626.

**Author details**

**6. References**

163-170.

Corresponding Author

 \*

D. Q. Truong and K. K. Ahn\*

Tokyo, Japan: 121.

	- [24] Kwok N.M, Ha Q.P, Nguyen T.H, Li J, Samali B (2006) A Novel Hysteretic Model for Magneto-Rheological Fluid Damper and Parameter Identification Using Particle Swarm Optimization. Sens. actuators A 132(2): 441-451.
	- [25] Kwok N.M, Ha Q.P, Nguyen M.T, Li J, Samali B (2007) Bouc-Wen Model Parameter Identification for a MR Fluid Damper Using Computationally Efficient GA. ISA Trans. 46(2): 167-179.
	- [26] Chang C.C, Roschke P (1998) Neural network modeling of a magnetorheological damper. J. intell. mater. syst. struct. 9(9): 755-764.
	- [27] Schurter K.C, Roschke P.N (2000) Fuzzy Modeling of A Magneto-Rheological Damper Using Anfis. Proc. of IEEE int. conf. on fuzzy sys.: 122-127.
	- [28] Wang D.H, Liao W.H (2001) Neural Network Modeling and Controllers for Magnetorheological Fluid Dampers. Proc. of 10th IEEE int. conf. on fuzzy sys.: 1323-1326.
	- [29] Wang D.H, Liao W.H (2005) Modeling and Control of Magneto-Rheological Fluid Damper Using Neural Networks. Smart mat. struc. 14(1): 111-126.
	- [30] Hanson B, Levesley M (2004) Self-sensing applications for electromagnetic actuators. Sens. actuators A. 116(2): 345-351.
	- [31] Mizuno T, Ishii T, Araki K (1998) Self-sensing magnetic suspension using hysteresis amplifiers. Con. eng. prac. 6(9): 1133-1140.
	- [32] Lama K.H, Chena Z.H, Nia Y.Q, Chanb H.L.W (2010) A magnetorheological damper capable of force and displacement sensing. Sens. actuators A. 158(1): 51-59.
	- [33] Ahn K.K, Truong D.Q, Thanh T.Q, Lee B.R (2008) Online Self Tuning Fuzzy PID Control for Hydraulic Load Simulator. Proc. of the IMechE - part I: j. sys. and con. eng. 222(2): 81-95.
	- [34] Truong D.Q, Ahn K.K (2011) Parallel Control for Electro-Hydraulic Load Simulator Using Online Self Tuning Fuzzy PID Technique. Asian j. con. 13(4): 522-541.
	- [35] Ahn K.K, Truong D.Q, Islam M.A (2009) Modeling of a Magneto-Rheological (MR) Fluid Damper using a Self Tuning Fuzzy Mechanism. J. of mech. sci. and tech. 23(5): 1485-1499.
	- [36] Truong D.Q, Ahn K.K (2010) Identification and application of black-box model for a self-sensing damping system using a magneto-rheological fluid damper. Sens. actuators A. 161(1-2): 305-321.
	- [37] Shames I.H, Cozzarelli F.A (1992) Elastic and Inelastic Stress Analysis. Prentice Hall, Englewood Cliffs, New Jersey 120-122 p.
	- [38] Gamota D.R, Filisko F.E (1991) Dynamic Mechanical Studies of Electrorheological Materials: Moderate Frequencies. J. of rheology 35(3): 399-425.


**4. Carbon sequestration in forests**

*Forest Degradation Around the World*

deterioration [14].

in the later succession.

grow [18].

**3. The role of natural regeneration**

of the land surface [6, 11]. Unfortunately, one-third of these original forests have been lost because of deforestation, resulting into large degraded areas. To date, the degraded areas are common in many countries, most of which are converted into other land uses, although some of them will regenerate into secondary forests [12]. Currently, the secondary forest can be found in all areas in the world and it is

The decrease of world's forest is caused by many activities such as clearing to open land for agriculture, roads, settlement, logging for timber, and cutting for fuel. Indirectly, forest is succumbing to the effect of environmental pollution and further threat by climatic change. Together, these causes have been responsible for decimating many of the world's forest, and they threaten to significantly degrade those that remain. Moreover, the burning of trees, shrubs, and other vegetation during land clearing and after logging further contributes to environmental

Sapling regeneration is important to restore forest formation after disturbance, balancing the loss of primary forest and as a main sequester of CO2. Their radial growth characteristics are the main performer of canopy closure, and it can be used to identify the canopy formation in the secondary forest after disturbances. Moreover, the canopy closure can be used to predict the vertical structure in the future; therefore, it will explain the pressure that will be given by dominant species

The continuing loss of primary tropical forest has made the secondary forest to be increasingly important for maintaining the ecology across the large forested landscapes [15]. The regeneration in the tropical secondary forest is strongly influenced by disturbances that resulted from natural tree fall in primary forest [15, 16]. Regeneration is widely recognized as being a key process in determining the structure and dynamism of tropical secondary forest [17]. For instance, the logging of mature trees can reduce seed resource in the forest [9]. Nevertheless, it may increase the light intensity, which reaches the forest floor and will stimulate the seeds' germination and growth of saplings. The ability of regeneration in the degraded areas is determined by the success of the seedlings and saplings to

The success of forest regeneration in terms of its distribution and diversity can be used as parameters to understand the level of forest succession and to figure out the population dynamic condition in the areas that have been controlled [18]. A lot of efforts have been made to study forest regeneration to provide information on ecological condition [18–20]; however, the data which have been collected cannot solve the same problem in different areas in the world. This is because every area in the world has different conditions which make it possible to create the different microclimate after fragmentation, and then influence the recovery process. The change of microclimate because of forest fragmentation will influence the plant species survival. Therefore, only species that can adapt to the environmental change will survive as different species have different capabilities to response to the ecological pressure. Furthermore, there are crucial chemical elements in the soil that have been recognized for growth of all plants [21], and the loss one of these elements will influence the plant growth, especially macronutrient. Besides the chemical elements, forest regeneration is influenced by light intensity which reaches under the canopy. The increasing of light intensity which reaches under the canopy will

estimated to reach about 64% of the total world forest [13].

**4**

Rapid development of industrial technologies in the world has contributed to the increment of greenhouse gases. The increase of greenhouse gases, such as CO2 in the atmospheres, has led to global warming phenomenon. Nowadays, rise of global temperatures has constituted a challenge to the mankind. The increment of CO2 today is estimated to be 1.3–1.9 ppm year−1 [22]. According to [23], from the analysis recorded in Malaysia, the average of temperatures has increased by 1.1°C within 50 years, which is consistent with the trend on global warming.

In reality, trees have important function in the carbon sequestration process [24–26]. In this case, the conservation of plant resources in a sustainable manner is the best solution. Forests play important roles in the CO2 sequestration process. Trees in the forests are important in order to perform the canopy closure in the forest and then influence the other species. Several studies have shown that the light intensity in the understory affects the growth and crown morphology of several tree species [17, 21, 27] and soil properties at allocation of plant biomass [28–31]. Additionally, the effect of soil properties and light intensity has been known to correlate with growth rate, and distribution pattern as well as species composition. In reality, each species has different capability to respond to their environmental factors.

In the tropical rainforests, nutrient cycling is probably one of the most significant ecological processes [32]. Similarly, according to Hamzah et al. [33], tropical rainforests are recognized as the richest ecosystems in the world in terms of structures, species diversity, and ability to maintain the global climate change by reducing the accumulation of greenhouse gases. According to Sands [34], tropical rainforests contain 70% of the world's plants and animals, 70% of the world's vascular plants, 30% of all birds, and more than 90% of all the invertebrates. There can be more than 200 tree species per hectare in tropical rainforest.


**1. Introduction**

Ocean waves contain large untapped renewable energy sources, which could be used to fulfil the energy demand [1]. Wave energy has the characteristics of high power density and wide distribution, and is renewable [2]. Hitherto, about one thousand wave energy converter (WEC) inventions have been patented [3], and more than 200 of them have entered the model testing stage [2]. However, present research mainly focuses on promoting the conversion efficiency of wave energy converters by improved mechanical design [4]. There are few research results to improve the wave energy generation efficiency from the perspective of increasing the wave energy density of the target sea area. Taking the characteristics of Bragg resonance into account is a good idea for wave energy focusing.

Bragg reflection or resonance is initially referred to as a special physical optical phenomenon of X-rays [5]. In the 1980s, Davies [6] first studied that Bragg resonance occurs if the bottom wavelength is an integral multiple of half wavelength of the incident waves. A followed idea originated for coastal protection was put forward. The strong reflection of water waves by seabed bars was investigated experimentally [7,8] as well as theoretically [9–11]. Afterwards, Bailard et al. [12] have extensively studied theoretically, computationally and experimentally to learn explicitly the feasibility of using longshore seabed bars to protect the shore and summarized its effectiveness and limitations. Some issues, including the partially standing-wave pattern [13] and shoreward increase

**Citation:** Zhang, H.; Tao, A.; Tu, J.; Su, J.; Xie, S. The Focusing Waves Induced by Bragg Resonance with V-Shaped Undulating Bottom. *J. Mar.Sci. Eng.* **2021**, *9*, 708. https:// doi.org/10.3390/jmse9070708

Academic Editors: Shih-Chun Hsiao, Wen-Son Chiang and Wei-Bo Chen

Received: 19 May 2021 Accepted: 24 June 2021 Published: 27 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

of wave amplitude [14–16] by Bragg reflection, were exposed. Howard & Yu [17] and Weidman et al. [18] have demonstrated that a specific phase relationship of wave-bottom interactions is significant, and sometimes it can cause exponentially varying standing-wave patterns. It is illustrated that Bragg resonance leads to offward intense reflection, even shoreward wave amplitude increase. Considering the strong effect of Bragg resonance, more studies focus on coastal protection or beach erosion from the aspect of the backward of the seabed bars.

Subsequent investigations of Bragg resonance mechanism in the field of coastal protection and engineering application pay more attention to the actual situations, e.g., oblique incident waves [10,19–22], sub- and superharmonic frequency [19–23], and seabed configuration [10,21,22,24].

Recently, Couston et al. [22] have proposed to revamp the Bragg resonance mechanism as a means of coastal protection by considering oblique seabed bars that divert, rather than reflect, shore-normal incident waves to the shore-parallel direction. Indeed, the novelty that with two superposed sets of oblique seabed bars reflecting waves along the shore-parallel direction, becomes efficiently deflected far to the sides, leaving a wake of decreased wave activity downstream of the patch. However, the wave characteristics of the forward of the oblique seabed bars are not studied by Couston et al. By contrast, few studies concern the wave reflection of the forward. From another perspective of the seabed bars reflecting waves, the wave-bottom interactions due to Bragg resonance may be of interest and be more in accordance with wave power generation need.

In the cases where closed-form solutions cannot be obtained, a wide range of numerical models are available to study water-wave scattering by seabed topographies [22]. These include the extended versions, e.g., the mild-slope equations [11,25], the coupledmode approach [26,27], the integral matching/discretized bottom method [28,29], the fully nonlinear Boussinesq equation adjusted for rapid bottom undulations [30], an asymptotic linear analytical solution (ALAS) [24], and the high-order spectral method [31]. Numerical investigations of Bragg scattering have helped explain several discrepancies between theory and experiments, including the difference between the observed and predicted Class II Bragg resonance frequency due to evanescent modes [32] and the resonant frequency downshift/upshift for the subharmonic/superharmonic Class III Bragg condition due to high-order nonlinearity [30].

Simultaneously, Bragg resonance theory has made some achievements in application fields. A series of submerged breakwaters have been developed to achieve coastal protection based on the Bragg resonance characteristic of reflecting waves to the sea in recent years. For the two-dimensional numerical cases, the parameters' influence on the Bragg resonance reflection coefficient has been investigated by scholars from different aspects, e.g., the shape, height, width, the number of the submerged breakwaters, the slope of the seabed, and wave conditions [33–37]. Generally speaking, these numerical simulations are based on various theoretical methods to study normally incident waves over a series of rectangular or trapezoidal breakwaters. Ning et al. [38] investigated numerically and experimentally the effect of the breakwater shape on hydrodynamic behavior in a two-dimensional flume. Shih & Weng [39] studied the interactions between waves and various combinations of undulating breakwaters in a three-dimensional basin. The reflection coefficient, transmission coefficient, and attenuation of wave energy were analyzed. Liu et al. [40] examined the Bragg reflection of waves by multiple submerged semi-circular breakwaters, by considering obliquely and normally incident waves independently. Although two-dimensional flume or three-dimensional basin experimental studies have been carried out, the study based on Bragg resonance for improving wave energy generation has not been widely considered. At present, research on the use of wave energy is gradually carried out. Water flow and wave propagation can be altered by artificial terrains or control devices. Zheng et al. [41] carried out a laboratory study on wave-induced setup and wave-driven current in a 2DH reef-lagoon-channel system. Elandt et al. [42] found that gravity waves can focus at a specific location by a concave mirror or a convex

lens of gravity waves. Tao et al. [43] investigated experimentally strong wave reflection induced by Bragg resonance. They found that wave reflection is effective to amplify the free-surface oscillation amplitude and focus the wave energy in front of the undulating bottom. Zhang & Ning [44] confirmed that reflected waves from the parabolic opening of the breakwater can travel towards a fixed focus position. For a specified wave environment, the wave heights at the focus positions can reach over several times of the incident wave heights, indicating wave energy multiplication. They all proposed the concept of the focal point. Instead of putting a multitude of wave energy harvesting devices over a large area, a large wave energy absorber can be placed at the focal point. However, the influence of peculiar spatial bottom layouts is ignored, which also affect the focal point and focusing area. Previous physical experiment research [43] shows that waves are strongly reflected in front of the undulating bottom. Based on that, this study considers the wave-focusing effect in a specific area by optimizing the undulating bottom pattern.

The objective of this study is to investigate both the characteristics of Bragg resonance induced by the regular waves and random waves. Three main research questions need to be answered: (1) How does the angle of inclination for the V-shaped undulating bottom affect the wave-focusing characteristics induced by the regular waves? (2) What analytical methods are used to determine the optimal angle of the V-shaped undulating bottom? (3) How do the steepness and spectrum width of random waves affect the wave-focusing characteristics with the V-shaped undulating bottom, and what is the law of their influence on wave height? By the high-order spectral (HOS) numerical simulation, these questions have been examined. The concept of V-shaped undulating bottom is first proposed, i.e., a horizontal V-shaped pattern with two continuous undulating bottoms inclined at the same angle perpendicular to the shoreline. The V-shaped undulating bottom can not only play the role of reflecting waves in front of the undulations, but also exploit the advantage of symmetrical "V-shaped", which can focus waves on the central axis and increase the wave heights in the focusing areas. However, waves in the ocean are stochastic, and the random waves contain wave components of different frequencies, each of which has different Bragg resonance effect due to interactions with the undulating bottom of fixed frequency. Most of the previous studies on wave Bragg resonance only consider the regular waves at the dominant resonance frequency. In addition, there are few studies on Bragg resonance involving the interactions between random waves and undulating bottom. In this work, we have further studied the Bragg resonance characteristics of the interactions between random waves and V-shaped undulating bottom by the numerical simulation of HOS. This study has been organized as follows. In Section 2, the generalized Bragg conditions, the high-order spectral (HOS) method, and the model establishment are introduced. In Section 3, Bragg resonance reflection coefficients of regular waves are studied from the perspectives of different angles *θ* and different values *f*. In Section 4, the Bragg resonance focusing characteristics of regular waves and random waves for V-shaped undulating bottom are expounded. In Section 5, the conclusions are provided.


*Article* **Identification of Aluminium Powder Properties for Modelling Free Air Explosions**

**Piotr W. Sielicki 1 , James Keith Clutter 2 , Wojciech Sumelka 1 , Tomasz Gajewski 1, \* , Michał Malendowski 1 , Piotr Peksa 1 and Robert Studzi ´nski 1**

**Abstract:** Aluminium is a component in many energetic formulations. Therefore, its combustion is one of the main thermochemical processes that govern the output from the energetics. Modelling aluminium combustion is a challenging task because the process is highly complex and difficult to measure. Here, tests of aluminium powder were conducted in an effort to isolate the burning of the aluminium and to determine an adequate representation of this process. Charges of 100 g and 500 g were tested, and the size of the Al/air cloud and the ratio of components in the Al/air mixture were determined, which has not been published previously. This information was used to assess the validity of the assumption that the detonation of the mixture was representative of the event. Parameters for the Jones–Wilkins–Lee equation of state for the explosive mixture and detonation products were defined. Simulations of the tests were performed, and the results were consistent with the field test data, indicating that detonation occurred when there was a mixture of 70–75% Al and 25–30% air by mass.

**Keywords:** aluminium powder; detonation; explosive; combustion; oxidation; equation of state

#### **1. Introduction**

The composition of an explosive is a key factor in its efficiency. Most standard military explosives and those used in research operations are condensed explosives, such as trinitrotoluene (TNT) and Plastic Explosive No. 4 (PE4). These substances are ideal for contact explosions due to their brisance and high-intensity shock wave propagation. Furthermore, they can be altered in different ways to change their explosive performances. One common alteration is the addition of various components, such as aluminium (Al) powder [1–3].

The addition of aluminium changes the energy release process and the overall output of an explosive at the microscopic level. Additionally, the ignition and combustion properties of an explosive can be modified, as shown by Liu et al. [4], by coating it with nanosized aluminium particles. Oxidation plays a crucial role in the combustion of aluminium. Gang et al. investigated the combustion and oxidation of Al nanoparticles at the atomistic level [5]. Depending on the amount of O2 and the temperature, oxidation can be grouped into three categories: mild, chain-like growth oxidation; moderate oxidation; and microexplosion-accelerated oxidation. The microexplosion-accelerated oxidation mechanism was investigated deeply by Gang et al. The burning of separated Al droplets in air was studied by Karasev et al. [6], as well as the mechanisms involved in alumina aggregate formation.

However, the microscopic level is not the only level of interest in terms of energy release processes; the combustion of aluminium nanoparticles has also been investigated at

**Citation:** Sielicki, P.W.; Clutter, J.K.; Sumelka, W.; Gajewski, T.; Malendowski, M.; Peksa, P.; Studzi ´nski, R. Identification of Aluminium Powder Properties for Modelling Free Air Explosions. *Materials* **2022**, *15*, 1294. https:// doi.org/10.3390/ma15041294

Academic Editor: Bolv Xiao

Received: 30 December 2021 Accepted: 6 February 2022 Published: 10 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the macroscopic level. For instance, Lewis et al. [7] studied hexahydro-1,3,5-trinitro-1,3,5 triazine (RDX) explosives with three types of aluminium nanoparticles. They showed that the addition of the nanoparticles can change the fireball temperature from 340 K to 4500 K. Gordon et al. [8] also studied the fireball of aluminised RDX, in addition to its shock wave energy. The authors reported that the shock energy is greater if aluminium is added to the high explosive rather than to the liner. Other studies of aluminised high explosives were conducted by Peuker et al. [9] and Carney et al. [10], in which optical methods were used. In the current study, the afterburn phase was analysed.

The contribution of Al combustion to the afterburn phase of aluminised explosives has been represented in various studies [11] using the Jones–Wilkins–Lee (JWL)–Miller model [12]. This model alters the energy term in the JWL equation such that it has the form

$$P = A \left( 1 - \frac{\omega}{R\_1 V} \right) e^{-R\_1 V} + B \left( 1 - \frac{\omega}{R\_2 V} \right) e^{-R\_2 V} + \frac{\omega (E + \lambda Q)}{V},\tag{1}$$

where *Q* is the additional heat released by the aluminium particle combustion, and *λ* is the progression variable indicating the degree to which the particle has reacted. The degree of reaction equation is given by

*dλ dt* = *a*(1 − *λ*) *m P n* , (2)

where *a* depends on the particle size, and *m* and *n* are reaction rate constants. This model represents the gradual addition of energy seen in the afterburn phase rather than the sudden release of energy observed in the aluminium powder explosives. Here, the rate of gasification, which is associated with the combustion rate equation, is relevant.

Other models have been introduced that also replicate the afterburn energy release due to Al combustion, but they assume that the reaction is proportional to the gasification of the particles [13]. This process is explicitly represented by the change in the diameter of the particles, which can be expressed as

$$\frac{d}{dt}\left[\frac{4}{3}\pi\left(\frac{D}{2}\right)^3\rho\right] = 4\pi\left(\frac{D}{2}\right)^2\dot{m}\_{\upsilon\nu} \tag{3}$$

which simplifies to

$$D = D\_0 - \frac{2\dot{m}\_v}{\rho}t,\tag{4}$$

where *D* is the particle diameter, *D*0 is the initial particle diameter, and *mv* is the gasification rate of the condensed material, which can be expressed as

$$
\dot{m}\_v = \frac{\sqrt{3}}{2\pi r^3} \sqrt{kT m} e^{-E\_v/RT} \,\text{,}\tag{5}
$$

where *r* is the radius of the molecule, *m* is the mass of each molecule, *k* is the Boltzmann constant (1.38 × 10−23 J/K), *T* is the temperature, *Ev* is the gasification enthalpy, and *R* is the gas constant (8.314 J/mol·K).

The afterburn combustion of Al can be used to redefine the degree of reaction, as follows

$$
\lambda = 1 - \left(1 - \frac{2\dot{m}\_v t}{\rho\_{Al} D\_0}\right)^3,\tag{6}
$$

where . *mv* is the gasification rate, *t* is time, *ρAl* is density of aluminum, and *D*0 is initial diameter.

Another group of models address the gasification of the Al particles using an empirical quasi-steady law [14]. The particle radius (*rp*) rate is

$$\frac{dr\_p}{dt} = -\frac{r\_p}{t\_b} \left( 1 + 0.276 \sqrt{Re} \right) \tag{7}$$

where *Re* is the Reynolds number, which is based on the relative velocity between the gas and the particle, and *tb* is the burning time based on

$$t\_b = K(D\_0)^2 \,\text{\,\,\,}\tag{8}$$

where *K* is an empirical constant. According to the literature, *K* is typically set equal to 4 × 106 s/m2 . The mass transfer from the solid to gaseous state is

$$\frac{dm\_p}{dt} = -\frac{d}{dt} \left(\frac{4}{3}\pi\rho\_p r\_p^3\right) \tag{9}$$

where the subscript *p* denotes the particle properties.

In the current study, the temperature of ignition depends on the state of the Al particles. If an oxide coating is present and is not cracked through some sort of physical process, then the melting point of Al2O3 would determine ignition. This melting point is approximately 2050 K. If the oxide layer is cracked, then the ignition temperature would be determined by the melting point of Al, which is approximately 950 K.

The temperature of the particles will depend on the heat feedback from the gas phase to the condensed phase. An investigation of aluminised explosives [15] has shown that the feedback can be represented using Fourier's law

$$q\_{\mathcal{E}} = \alpha\_{\mathcal{S}} \left( \frac{1}{r} \right) \left( \frac{dT}{dt} \right)\_{s'} \tag{10}$$

where *αg* is the thermal conductivity, *r* is the burning rate, and *s* denotes the gradient at the surface.

Hence, an understanding of the energy released during aluminium combustion is key to optimising the design of energetic systems and would allow the quantification of aluminium's effect on structures. In this paper, aluminium powder was analysed at the macroscopic level; this approach enables the straightforward development of online (fast) tools that can predict the behaviour of aluminium powder explosives. The current objective was to develop an adequate representation of Al combustion to provide loading predictions for use in the analysis of the interaction between munitions and structures. From this point of view, the robust modelling of such explosives is no trivial task.

Equations of state (EOS) are a key concept in modelling this class of energetics. Various EOS approaches are found in the literature. They have been proposed for aluminised explosive products [16], along with more popular approaches, such as the JWL formula [17–19]. The JWL approach is well known and has the broadest application; thus, it was adopted in the current study. Here, the expansion process of the detonation products from the composite energetics was determined to be crucial. This expansion process affects the prediction of the pressure loading at a particular distance, which is an essential characteristic in determining the effect of the explosive.

This paper analysed the combustion of Al, as this critical process affects the energy outputs of aluminised explosives. The oxidiser that was used in the tests varied according to the specifics of the energetics. In most explosives, the oxidiser is obtained from the products of the organic reaction that forms CO, CO2, and H2O. It has been shown that Al combustion is determined by the concentrations of CO2 and H2O [20]. Therefore, Al combustion often occurs before any ambient oxidiser is introduced to the Al. Here, pure Al was used to ensure that the oxidation process was related to the mixing of Al powder with ambient air. This provided insights into the combustion of Al, even when the oxidiser originated from typical detonation products.

In this study, a series of range detonation tests with pure aluminium powder mixtures were performed, which has not been published previously. The mass of the charges varied from 0.1 to 0.5 kg. The primary goal of the tests was to measure the pressure–time histories, and these results were replicated using modelling tools. Seven tests were performed, and

full data recordings were obtained at standoff distances of approximately 1 and 2 m. The tests were replicated using a modelling tool that incorporated the effects of Al combustion. The EOS and reaction parameters were first defined using thermodynamic codes and analytical tools. Subsequently, as described in the final part of this paper, these values were altered to match the test results.

The paper is structured as follows. In Section 2, the test details, computational assumptions, and an explosive governing equation are presented. Section 3 identifies the explosive parameters for the 100 g charge and presents their validation for the 500 g charge. Section 4 presents the conclusions of the study.

#### **2. Materials and Methods**

#### *2.1. Field Test Measurements*

In this section, data from the field tests for the 100 g and 500 g aluminium powder charge detonation are described. The experiments were conducted so as to represent a free air blast. Data from the 100 g detonation were used to determine the aluminium modelling parameters. The tests were performed using 6 µm of Al powder with an estimated density of 2.2 to 2.7 g/cc. The aluminium powder is shown in Figure 1a. Photographs from the test site are presented in Figure 1b,c, including a snapshot of the explosion (Figure 1b) and an overview of the explosion area (framed by wooden poles), along with the high-speed camera (Figure 1c). Figure 1d presents details of the explosive devices' mounting scheme. The aluminium powder was hung in a foil bag attached to a string and was located 100 cm from the ground (L) directly between two wooden poles that were 200 cm high (H). The distance between the poles was 500 cm (W). Blast pressure pencil probes were placed at heights of 100 cm (b) and 200 cm (a) from the charge.

Seven detonations were conducted during the field tests. In the first two detonations, the aluminium powder charge mass was 100 g; in the next two, the charge mass was 200 g, and in the last three, the charge mass was 500 g. The data from the second 200 g detonation were rejected. In this case, during aluminium combustion, the foil bag was damaged, and the powder poured out of the bag. Thus, the data could not be considered for further analysis. Because the 200 g case was not represented by two detonations, this charge mass was not taken into account in the subsequent analyses.

**Figure 1.** *Cont*.

(**c**)

**Figure 1.** Test set up: (**a**) example of tested mixture of aluminium powder; (**b**) snapshot of the explosion; (**c**) overview of the explosion area; and (**d**) explosive devices' mounting scheme.

The pressure–time histories from all tests were measured at 1 and 2 m from the explosive simultaneously. Two ICP® blast pressure pencil probes were used simultaneously to acquire the data—these are the same probes that were used in [21]. The maximum pressure limit for these sensors is 345 kPa. The results are shown in Figures 2 and 3. All measurements are presented, including those that were rejected (i.e., pressure histories no. 3 (200 g) and no. 4 (labelled as Err, i.e., error)).

− The detonations were filmed with a high-speed camera to obtain deeper insights into the reaction process to ensure adequate modelling assumptions. In these types of experiments, a high-speed camera is typically used to assess the deformation due to the blast wave [22] and determine the speed of projectiles and fragments, e.g., [23,24]. In our tests, the high-speed camera was used to estimate the volume of the mixed powder. Figure 4 shows sample images captured immediately after detonation. The time after ignition is shown next to the sequence of images. Powder burning at the time of ignition and thereafter can be observed.

57

**Figure 2.** Overpressure measurements at 1 m.

**Figure 3.** Overpressure measurements at 2 m.

**Figure 4.** Sequence of images from the 100 g test. The time presented refers to the time after ignition.

4 + 3 ⟶ 2ଶଷ

If a spherical charge is assumed and the Al density is 2.7 g/cc, the volume of the initial charge would be 37.04 cc, and the charge would have a diameter of 2.07 cm. Notably, there was an apparent spike in output at the 18.1 ms mark. The Al/air cloud at the 18.1 ms mark had a diameter of approximately 41.4 cm and a volume of 37.15 cc. This resulted in an Al density of 2.69 × 10−3 g/cc. Taking the air to be at atmospheric pressure, the volume of powder and air would contain a mixture of 56% Al and 44% air considering masses. This information was used to derive an adequate representation of the explosion process.

The overall pattern for the 500 g tests was similar to that of the 100 g tests. For the 500 g charge, the "explosion" (rapid release of energy) appeared to occur between 2.3 and 5.8 ms after ignition. The diameter of the Al/air cloud at that time varied from 40 to 54 cm. Figure 5 shows an image taken after the "explosion" of the 500 g Al charge. More burning particles of Al were expelled outward after the "explosion" (white arrows) of the 500 g charge compared to that of the 100 g charge, suggesting that the 500 g charge was less effective than the 100 g charge. The analysis of the 500 g test data allowed an adequate representation of the Al combustion to be determined.

**Figure 5.** Image of the 500 g test (Test 7), in which burning Al particles (white arrows) were expelled.

#### *2.2. Numerical Modelling Assumptions*

The computational model that was developed to simulate the energetics used a Cartesian Adaptive Mesh (CAM) framework; this model simulates scenarios involving a blast, explosion, and release of materials (thus, CAMBER). CAMBER is an object-oriented framework that utilises a variety of material models and reaction laws. The mesh adaptation was used to resolve structures under severe deformation, and the code maps out the locations of gradients in properties [25]. In CAMBER, a finite-volume formulation is used. Various flux calculation schemes have been assessed and proven to be useful in multi-material flows [26]. The advection upstream splitting method (AUSM+) scheme was selected for its simplicity and robustness [27]; high spatial order is achieved via a monotonic upstreamcentred scheme for conservation law (MUSCL) extrapolation with limiters. Given the small time steps required to resolve the evolving flow field, explicit time integration was used. The code can be applied in 2D, 2D-axisymmetric, or 3D modes.

#### *2.3. Governing Equation for Explosions*

The objective of this study was to formulate an approach to simulate explosions involving aluminium. The energetics used in the present research were for pure aluminium, so the reaction causing the energy release was as follows

4 + 3 ⟶ 2ଶଷ 4Al + 3O → 2Al2O3. (11)

Hence, if Al is part of the explosive composition, the oxidation could arise from the CO2 and H2O in the detonation product. In this case, the oxidiser was the O2 in the air. There was clearly burning Al prior to the "explosion" at 18.1 ms (see Figure 4), but at this time point, there was a sudden release of energy. It can be assumed that until this time point, the percentage of Al burned was negligible. Furthermore, the assumption was made that the explosive cloud was 56% Al and 44% air by mass (as described in Section 2.1).

A thermochemical analysis was performed as a first estimate; it was assumed that the reaction was essentially the detonation. The Chapman–Jouguet (CJ) conditions are listed in Table 1. A JWL EOS of the following form

$$P = A\left(1 - \frac{\omega}{R\_1 V}\right)e^{-R\_1 V} + B\left(1 - \frac{\omega}{R\_2 V}\right)e^{-R\_2 V} + \frac{\omega E}{V},\tag{12}$$

was used for the explosive material and detonation products. The parameters used are listed in Table 2. The JWL parameters were defined for the Al/air mixture based on the behaviour of such mixtures under shock loading. The rate of the reaction was derived from the detonation velocity.

**Table 1.** Detonation conditions assumed in the computational thermochemical code for the Al/air explosion.

**Table 2.** JWL parameters assumed in the computational thermochemical code for the detonation products.

\* This value was adjusted to 2.7 × 10−3 after initial simulations of the 100 g case to better match the test data.

#### **3. Identification and Validation**

#### *3.1. Identification—Simulation of 100 g Tests*

The parameters defined in Section 2 were used to simulate the 100 g tests and identify a set of explosive parameters. The initial condition was that there was a 37.04 cc Al/air cloud. Figure 6 shows the sequence of images from the simulation. The contour maps for the concentration of reacting mixture (Al/air), the reaction products, the density, and the pressure are shown at four selected time points. The time noted is relative to *t* = 0 for the tests. Figure 7 shows a comparison between the test data for the pressure–time histories recorded during the simulation at the 1 and 2 m locations.

**Figure 6.** Sequence of images from the 100 g test simulation: (**a**) Al/air mixture concentration; (**b**) detonation product concentration; (**c**) density; and (**d**) pressure are shown at four time points.

The initial pulses in pressure at both locations were similar to the test data in both magnitude and duration. The secondary pulses due to ground reflection were somewhat lower in magnitude and were delayed compared to the test data. This difference could be due to slight differences in the heights. There was an intermediate pulse in one of the 1 m test recordings. While the cause of this pulse was not clear, it was not seen in the simulation results.

These results suggest that the predominant release of energy from the combustion of the Al/air cloud can be represented as a sudden event, or, as we label it here, the "explosion". It appears that the combustion event was essentially a mixing-controlled process. The explosion did not occur until the mass of Al powder expanded to the point where there was sufficient oxidation. Based on the size of the cloud at the time of the explosion and the mass of the Al powder, the cloud distribution was approximately 70% Al and 30% air by mass, providing an oxidiser-to-fuel mole ratio (moles O/moles Al) of 0.17, far below the stoichiometric ratio for the reaction in Equation (1).

**Figure 7.** Comparisons between the test data and simulated values for the pressure versus time history at the 1- and 2 m locations for the 100 g charge.

#### *3.2. Validation—Simulation of the 500 g Tests*

The 500 g tests were simulated in the validation stage. The same modelling parameters were used in terms of EOS and the rate of reaction. First, the test imagery was used to estimate the size of the spherical Al/air cloud at the point at which there was a noticeable sudden release of energy (i.e., the "explosion"). Three videos were reviewed, and the time of the explosion varied from approximately 3 to 6 ms. The size of the cloud was estimated to be approximately 52 cm in diameter. This diameter would result in a cloud composed of 84% Al and 16% air by mass. Although the fuel-to-oxidiser ratio changed slightly from the 100 g case, the same parameters used in that case were applied. The pressure–time histories predicted at the 1 and 2 m locations are shown in Figure 8. At both locations, the peak simulated pressures were lower than the peak measured pressures.

**Figure 8.** Comparison between test data and simulated values for the pressure versus time history for the 500 g Al/air cloud with a 52 cm diameter.

Because the cloud estimated from the video produced a lower output, several cloud sizes were applied until a good match between the measured and simulated pressure–time histories was achieved. Figure 9 shows the results when the diameter of the Al/air cloud was assumed to be 63 cm. The comparison is moderately good, wherein the peaks at the 2 m location are slightly lower than those in the test data. Using the assumed size of the cloud and the mass of Al, the cloud was found to be 75% Al and 25% air by mass. This result was similar to the ratio found for the 100 g results. The study presented shows that combining the experimental and advanced numerical approaches creates the ability to obtain the synergic effect for better understanding the physical phenomena, similar to [28].

**Figure 9.** Comparison between test data and simulated values for the pressure versus time history for the 500 g Al/air cloud with a 63 cm diameter.

#### **4. Conclusions**

Combustion is the key process in many energetic systems. In this paper, the combustion of Al was investigated using a series of tests and simulations, which has not been published previously in the literature. Here, the Al was isolated and allowed to react with just oxygen in ambient air; the oxidation arose from the detonation products. The results of the 100 g and 500 g Al charges indicated that the sudden explosions occurred when a mixture of 70–75% Al and 30–25% air, by mass, was obtained.

Among others, the novelty of the paper is that the output from the explosions was replicated using the JWL state equation for the product and assuming detonation of the mixture. The detonation was modelled using a prescriptive method that set the burn rate based on the detonation velocity. The delay from the initiation of the event until the time of the explosion is currently under investigation. Further studies will aim to determine whether this delay is related to the time required for the gasification of the Al or the event is a mixing-controlled process.

**Author Contributions:** Conceptualisation, J.K.C. and P.W.S.; methodology, J.K.C. and P.W.S.; software, P.W.S.; validation, J.K.C., R.S., W.S. and P.W.S.; formal analysis, J.K.C., R.S., W.S., T.G. and P.W.S.; investigation, J.K.C. and P.W.S.; resources, M.M. and T.G.; data curation, M.M. and T.G.; writing—original draft preparation, J.K.C., P.W.S. and T.G.; writing—review and editing, P.W.S., R.S., T.G., M.M., W.S. and P.P.; visualisation, T.G., M.M. and P.P.; supervision, W.S., T.G. and P.W.S.; project administration, P.W.S. and T.G.; funding acquisition, P.W.S., W.S. and P.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Centre for Research and Development Poland under the grant DOB-BIO10/01/02/2019 within the framework of the Defence and Security Programme.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

#### **References**