**Preface to "Characterization and Processing of Complex Materials"**

This reprint provides a forum for scientists and engineers to share and discuss their pioneering/original findings or insightful reviews on the "Characterization and Processing of Complex Materials". Reports on successfully coupling complex industrial materials/waste/natural ores with non-destructive characterization method(s) and the unique processing developments were particularly encouraged. On the other hand, any kind of contribution under the broad framework of the understanding and production of complex materials was also highly regarded.

This reprint gathers papers on the characterization and processing of complex materials with the following keywords: non-destructive characterization; analytical chemistry and instrumentation; waste management and recycling; selective material recovery/processing; sustainability and circular economy. This includes the waste management and recycling of complex materials (e.g., food and electronic waste) through the development of a new characterization method and/or process intensification. For example, the proper characterization of complex/heterogeneous materials is still an extremely challenging task, since the majority of characterization methods analyze either the average characteristics of the whole material or a narrow area of specific interest. On the other hand, the efficient processing/recycling of complex materials/waste should be strongly supported by their characterization from the viewpoint of economic and environmental concerns.

> **Akira Otsuki** *Editor*

### *Editorial* **Editorial for the Special Issue: "Characterization and Processing of Complex Materials"**

**Akira Otsuki 1,2,3**


#### **1. Introduction**

The Special Issue aimed to provide a forum for scientists and engineers to share and discuss their pioneering/original findings or insightful reviews on the "Characterization and Processing of Complex Materials". Reports on the achievement of coupling the complexity of industrial materials/wastes/natural ores with non-destructive characterization method(s) and the unique advancement of their processing were particularly encouraged. On the other hand, any kind of contribution under the broad framework of the understanding and production of complex materials was also highly regarded.

The Special Issue was intended to gather papers on the characterization and processing of complex materials with the following keywords: non-destructive characterization; analytical chemistry and instrumentation; waste management and recycling; selective material recovery/processing; sustainability and circular economy. One of them is the waste management and recycling of complex materials (e.g., foods and electronic wastes) via the development of a new characterization method and/or process intensification. For example, the proper characterization of complex/heterogeneous materials is still an extremely challenging task, since the majority of characterization methods analyze either the average characteristics of the whole material [1] or a narrow area of specific interest [2,3]. On the other hand, the efficient processing/recycling of complex materials/wastes should be strongly supported by their characterization from the viewpoint of economic and environmental concerns [2].

#### **2. The Special Issue**

The Special Issue contains 10 papers from around the globe on the topics related to the "Characterization and Processing of Complex Materials", and specifically discusses the following contents.

Wu et al. [4] and Zhou et al. [5] reported their nanobubble studies. Wu et al. [4] developed a method to control the crystalline and morphology of calcium carbonate with the aid of nanobubbles and showed the structure transition from vaterite to calcite. They explained that nanobubbles can coagulate with calcite, work as a binder between calcite and vaterite and accelerate the transformation from vaterite to calcite. They also explained their proposed mechanism by the potential energy calculated using the extended DLVO (Derjaguin–Landau–Verwey–Overbeek) theory. Zhou et al. [5] studied five kinds of nanobubbles (N2, O2, Ar + 8%H2, air and CO2) in deionized water and a salt aqueous solution prepared using the hydrodynamic cavitation method. They measured the mean size and zeta potential of the nanobubbles using a light scattering system, while the pH and Eh of the nanobubble suspensions were measured as a function of time. The nanobubble stability was predicted and discussed by the total potential energies between two bubbles by the extended Derjaguin–Landau–Verwey–Overbeek (DLVO) theory.

**Citation:** Otsuki, A. Editorial for the Special Issue: "Characterization and Processing of Complex Materials". *Materials* **2023**, *16*, 3830. https:// doi.org/10.3390/ma16103830

Received: 18 May 2023 Accepted: 18 May 2023 Published: 19 May 2023

**Copyright:** © 2023 by the author. 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/).

Fukushima et al. [6,7] reported their studies on the modeling and simulation of electrical disintegration. Fukushima et al. [6] developed a simulation method using equivalent circuits of granite to better understand the crushing process with high-voltage (HV) electrical pulses that have been previously identified as a selective liberation method. Using their simulation works, they calculated the electric field distributions, produced heat and temperature changes in granite when an electrical pulse was applied. Their results suggested that the void volume in each mineral is important in calculating the electrical properties and confirmed that considering both the electrical conductivity and dielectric constant of granite can more accurately represent the electrical properties of granite under HV electric pulse application. Fukushima et al. [7] prepared equivalent circuit models in order to represent the electric and dielectric properties of minerals and voids in a granite rock sample, considering HV electrical pulse application. Using the electric and dielectric properties and the mineral compositions measured, they prepared ten patterns of equivalent circuit models. The values of the electric circuit parameters were determined from the known electric and dielectric parameters of the minerals in granite. The average calculated data of the electric properties of granite agreed with the measured data.

Elphick et al. [8] applied their scanning electron microscopy (SEM)-based non-destructive imaging technique for the characterization of nanoparticles synthesized using X-ray radiolysis and the sol-gel method. They observed the interfacial conditions between the nanoparticles and the substrates by subtracting images taken by SEM at controlled electron acceleration voltages to allow backscattered electrons to be generated predominantly below and above the interfaces. They found the interfacial adhesion to be dependent on the solution pH used for the particle synthesis or particle suspension preparation, proving the change in the particle formation/deposition processes with pH as anticipated and in agreement with the prediction based on the DLVO theory.

Kanari et al. [9–12] reported their characterization results of chromite concentrate [9], copper anode slime [10,12] and alkali ferrates synthesized from industrial ferrous sulfate [11]. Their characterization consists of inductively coupled plasma atomic emission spectroscopy (ICP-AES), SEM-EDS and X-ray diffraction in order to extract information on the chemical composition, morphology and mineralogical composition of their products prepared under different conditions.

Zhong and Xu [13] reported high-temperature mechanical behaviors of SiO2-based ceramic cores for the directional solidification of turbine blades. They conducted isothermal uniaxial compression tests of ceramic core samples on a Gleeble-1500D mechanical simulator with an innovative auxiliary thermal system and investigated the stress–strain results and macro- and micro-structures of SiO2-based ceramic cores. They characterized the microstructures using SEM. Based on their experimental data, they established a nonlinear constitutive model for high-temperature compressive damage. The statistical results of the Weibull modulus showed that the stability of hot deformation increases with the increase in temperature. They confirmed the applicability of their model by comparing their results between the predictions of the nonlinear model and the experimental values.

**Acknowledgments:** The guest editor wishes to thank all the authors, reviewers and the editors who contributed to this Special Issue.

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

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

### *Article* **Long-Term Stability of Different Kinds of Gas Nanobubbles in Deionized and Salt Water**

**Yali Zhou 1, Zhenyao Han 1, Chunlin He 1, Qin Feng 1, Kaituo Wang 1, Youbin Wang 1, Nengneng Luo 1, Gjergj Dodbiba 2, Yuezhou Wei 1, Akira Otsuki 3,4 and Toyohisa Fujita 1,\***


**Abstract:** Nanobubbles have many potential applications depending on their types. The long-term stability of different gas nanobubbles is necessary to be studied considering their applications. In the present study, five kinds of nanobubbles (N2, O2, Ar + 8%H2, air and CO2) in deionized water and a salt aqueous solution were prepared by the hydrodynamic cavitation method. The mean size and zeta potential of the nanobubbles were measured by a light scattering system, while the pH and Eh of the nanobubble suspensions were measured as a function of time. The nanobubble stability was predicted and discussed by the total potential energies between two bubbles by the extended Derjaguin–Landau–Verwey–Overbeek (DLVO) theory. The nanobubbles, except CO2, in deionized water showed a long-term stability for 60 days, while they were not stable in the 1 mM (milli mol/L) salt aqueous solution. During the 60 days, the bubble size gradually increased and decreased in deionized water. This size change was discussed by the Ostwald ripening effect coupled with the bubble interaction evaluated by the extended DLVO theory. On the other hand, CO2 nanobubbles in deionized water were not stable and disappeared after 5 days, while the CO2 nanobubbles in 1 mM of NaCl and CaCl2 aqueous solution became stable for 2 weeks. The floating and disappearing phenomena of nanobubbles were estimated and discussed by calculating the relationship between the terminal velocity of the floating bubble and bubble size.

**Keywords:** extended DLVO theory; mean size; zeta potential; Ostwald ripening effect; Stokes equation

#### **1. Introduction**

Nanobubbles have some unique properties, unlike conventional milli- to micro-bubbles, such as high mass transfer [1], long-term stability [2–10], high zeta potential, high surface to volume ratio, and generating free radicals when collapsing [11,12]. Nanobubbles can be divided into surface nanobubbles absorbed on solid surfaces and bulk nanobubbles dispersed in aqueous solutions, experiencing Brownian motion. Bulk nanobubbles have diameters of less than 1 micrometer [13]. Because of their unique physico-chemical properties, nanobubbles can be used in various application fields, e.g., improvement of plant growth and productivity [14], membranes cleaning [15–17], waste-water treatment [1,18–22], visualization improvement as the ultrasound contrast agent [23], froth flotation [24,25], improvements of methane production in the anaerobic digestion [26,27], applications in food processing [28,29] and reactions with concrete using CO2 nanobubbles [30,31].

**Citation:** Zhou, Y.; Han, Z.; He, C.; Feng, Q.; Wang, K.; Wang, Y.; Luo, N.; Dodbiba, G.; Wei, Y.; Otsuki, A.; et al. Long-Term Stability of Different Kinds of Gas Nanobubbles in Deionized and Salt Water. *Materials* **2021**, *14*, 1808. https://doi.org/ 10.3390/ma14071808

Academic Editor: Daniela Kovacheva

Received: 8 March 2021 Accepted: 2 April 2021 Published: 6 April 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/).

Different kinds of nanobubbles have different application potentials. H2 nanobubble gasoline blends can improve combustion performance, compared with conventional gasoline [32]. A N2 nanobubble water addition can enhance the hydrolysis of waste activated sludge and improve methane production in the process of anaerobic digestion [26]. O2 nanobubbles produce the methane in anaerobic digestion of cellulose [33] and CO2 bulk nanobubbles can be used in food processing [29].

Our group also reported free radical degradation in water using different kinds of nanobubbles, i.e., H2 in Ar, O2, N2, CO2 and a mixture of H2 in Ar and CO2. The hydroxyl radical was scavenged, and the superoxide anion was diminished by mixing the carbon dioxide nanobubbles after hydrogen nanobubbles existence in the water [11]. The antioxidant effect of H2 nanobubbles in water was found to suppress tumor cell growth [12]. In those applications, it is important to investigate the stability and quality of nanobubbles in water, depending on time for further utilization.

Since there are so many potential applications, it is important to study fundamental characteristics and properties of bulk nanobubbles in complex solutions. The low concentration of 1 mM (milli mol/L) of NaCl could stabilize O2 micro- and nano-bubbles for at least one week [7] because there was a shielding effect in low concentrations of NaCl. On the other hand, the higher concentration of NaCl decreased the nanobubble concentration more quickly [7]. Nirmalkar et al. (2018) performed the calculation of the interaction energies between air nanobubbles in 0–20 mM NaCl solutions by using the Derjaguin–Landau–Verwey–Overbeek (DLVO) theory. According to their calculation, the interaction potential energy was positive above pH 4 in deionized water and it decreased with an increasing salt concentration. Beyond a certain critical concentration of NaCl (between 10 to 20 mM), the system becomes unstable [8]. Meegoda et al. (2019) calculated the electrical double layer potential between O2 nanobubbles in 0–0.1 M NaCl solutions and reported the stability of nanobubbles in 1 mM NaCl concentration [9]. Hewage et al. (2021) also reported the stability of air nanobubbles in 1 mM of different ion valence electrolyte solutions (i.e., NaCl, CaCl2, FeCl3) was confirmed over one week [10]. However, these research studies are only partial. In order to have a better and fuller understanding of nanobubble characteristics and have a reasonable interpretation of nanobubble behaviors in different solution environments, it is important to have a comprehensive study.

In the current study, five different kinds of nanobubbles of N2, O2, 8% H2 in Ar, CO2 and air were compared with their existence period, average size, and zeta potential, as well as their suspension pH and Eh as a function of time to be considered in several applications. The property changes of nanobubbles in the presence of salt as a function of time were observed and discussed. The stability of nanobubbles via our experimental studies was interpreted and discussed by using the extended DLVO theory.

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

#### *2.1. Materials*

Deionized water with a resistivity of 18.2 MΩcm prepared by the Classic Water Purification System from Hitech instruments CO., Ltd. (Shanghai, China) was used for all the experiments. Different gases (N2, O2, Ar + 8%H2 (8% H2 and 92% Ar mixed gas), air and CO2) were used to prepare nanobubble water. Sodium chloride (NaCl), calcium chloride (CaCl2), aluminum chloride hexahydrate (AlCl3.6H2O) were used to prepare 1 mM salt solutions and nanobubble water in order to study the nanobubbles' stability in the case of salt addition, considering the Schultze Hardy rule, i.e., *CCC* ∝ <sup>1</sup> *<sup>z</sup>*<sup>6</sup> where *CCC* is the critical coagulation concentration, and *z* is the ionic valence [34]. Dilute HCl or NaOH was used to adjust the pH of the nanobubble water. A 10 L glass container was used to store the nanobubble water temporarily and 50 mL plastic centrifuge tubes were used to store the nanobubble water after they were rinsed well with deionized water.

#### *2.2. Generation of Different Kinds of Nanobubbles*

Nanobubble water was produced by a ultra-micro bubble generator (XZCP-K-1.1) (100 nm to 10 μm) provided by Xiazhichun Co. Ltd. (Kunming, China). When the generator was turned on, gas from a gas cylinder was pumped into the machine in a negative pressure, and inlet and outlet pipes were immersed in deionized water/electrolyte solutions in the 10 L glass container where there was a circulation between the inlet and outlet, as shown in Figure 1A. The machine mixed the gas and liquid, and then high-density, uniform and "milky" nanobubble water (Figure 1B) was produced through a nozzle by the hydrodynamic cavitation method. The cloudy and milky nanobubble water gradually became clear through the process of the microbubbles rising and collapsing at the air–water interface (Figure 1C). This process took 2 to 3 min (Figure 1D). The machine worked for 15 min to generate high number density nanobubble water of more than 10<sup>8</sup> bubbles/mL (obtained from NanoSight, NS300, Malvern (Worcestershire, UK) outsourcing).

**Figure 1.** Procedure of nanobubble generation. (**A**) Before generation; (**B**) During generation of micro-nano bubbles; (**C**) Stop the generation of bubbles; (**D**) After standing for 2 to 3 min.

After the above-mentioned preparation steps, the nanobubble water had a temperature of about 313 K. It was then left to cool down to room temperature. Next, from the prepared nanobubble water, the large sizes of the bubbles, from 1 to 10 μm, were removed by centrifugal treatment in the following manner. The nanobubble water was stored in a 50 mL centrifuge tube; then, centrifugal treatment was performed in 6000 rpm, i.e., 31.4 m/s peripheral velocity, for 6 min to remove possible impurities and large bubbles. After centrifugation, the size distribution, zeta potential, Eh and pH were measured from one centrifuge tube to record the first day's data while other sealed centrifuge tubes with nanobubble water were kept at room temperature (298 K) for continuous measurements.

#### *2.3. Characterization of Bulk Nanobubble Suspensions*

The nanobubble size was measured by the dynamic light scattering method (DLS, NanoBrook Omni, Brookhaven Instruments, Holtsville, NY, USA). The machine's measurement range is between 1 nm and 10 μm. Since nanobubbles experience Brownian motion, the scattered light fluctuates as a function of time due to the particles' random movements, and the diffusion coefficient can be obtained from a computer digital correlator. The particle hydrodynamic diameter can be calculated by the Stokes–Einstein equation [35] as described below:

$$D\_T = \frac{KT}{3\pi\eta d\_h} \tag{1}$$

where *DT*, *k*, *T*, *η*, *dh* is the diffusion coefficient, Boltzmann constant, liquid absolute temperature, viscosity, and hydrodynamic diameter, respectively.

One measurement duration time was 3 min, and the instrument performed the size measurement 3 times. After size measurements, the zeta potential was measured by the same instrument that performed the 3 measurements, consisting of 20 runs/measurement. Every experiment was performed in duplicate and then the average size and measurement error were calculated and plotted.

The zeta potential value was obtained through the micro-electrophoresis method by the phase analysis light scattering method (NanoBrook Omni, Brookhaven Instruments). The zeta potential of the bubbles in the aqueous solution with salt was calculated by using the Smoluckowski model *f*(*κa*) = 1 in Equation (2) for *a* >> 1/*κ* (*κa* >> 1); however, the zeta potential of the bubbles in deionized water was calculated by using the Hückel model *f*(*κa*) = 2/3 in Equation (2) individually for *a* << 1/*κ* (*κa* << 1) [34] because of the small bubble size. The electrophoretic mobility *μ* in the Henry equation is shown as follows:

$$
\mu = \frac{\varepsilon\_r \varepsilon\_0 \tilde{\zeta}}{\eta} f(\kappa a) \tag{2}
$$

where *εr*, *ε*<sup>0</sup> are relative permittivity and permittivity of free space, respectively, *ζ* is the zeta potential, *η* is fluid viscosity, *κ* is the Debye length, *a* is bubble radius and *f*(*κa*) is the Henry function [36].

The viscosity, refractive index and relative dielectric constant of water at 293 K are 1.002 × <sup>10</sup>−<sup>2</sup> poise, 1.332, and 80.2, respectively [37]. Nanobubble water pH and Eh values were measured by pH meter (Inlab Expert Pro, Mettler Toledo, Greifensee, Switzerland) and Eh meter (by oxidation-reduction potential meter, Inlab Redox, Mettler Toledo), respectively.

#### **3. Results and Discussion**

#### *3.1. Zeta Potential of Different Kinds of Nanobubbles*

The zeta potentials of prepared different kinds of nanobubbles (N2, O2, Ar + 8%H2, air and CO2) depending on pH were measured and are plotted in Figure 2. The pH was adjusted from the natural pH by adding NaOH or HCl to achieve the desired pH. The natural pH of nanobubbles in deionized water was between a pH of 6 to 7, except for CO2 nanobubble suspensions. The zeta potential of N2, O2, Ar + 8%H2 nanobubbles were negative at about −15 to −30 mV at a natural pH of 6 to 7, while the zeta potential of CO2 nanobubbles in deionized water was positive at about +10 mV at a natural pH of 4 to 4.5. The solubility of gas, dielectric constant of gas [37] and the isoelectric point (IEP) of nanobubbles determined in this study are listed in Table 1. Since the solubility of CO2 is high (7.07 × <sup>10</sup>−4, Table 1), the HCO3 − adsorption on the CO2 bubble surface can cause a positive zeta potential in the deionized water due to the concentration of counter ions. The CO2 solubility phenomena can be explained in the following Equations (3)–(8) [38] assuming that CO2 gas is dissolved in water.

$$\rm CO\_2 + H\_2O \rightleftharpoons H\_2CO\_3 \tag{3}$$

$$\frac{H\_2\text{CO}\_3}{\text{CO}\_2} = 1.7 \times 10^{-3} \text{ at } 298 \text{ K} \tag{4}$$

$$H\_2\text{CO}\_3 \rightleftharpoons H^+ + H\text{CO}\_3^- \tag{5}$$

$$K\_{a\_1} = \frac{[H^+] \left[HCO\_3^-\right]}{[CO\_2] + [H\_2CO\_3]}, \text{ } pK\_{a\_1} = 6.35\tag{6}$$

$$\rm{HCO\_3^-} \rightleftharpoons \rm{H^+} + \rm{CO\_3^{2-}} \tag{7}$$

$$K\_{a2} = \frac{[H^+]\left[CO\_3^{2-}\right]}{\left[HCO\_3^-\right]}, \ pK\_{a2} = 10.33\tag{8}$$

**Figure 2.** Zeta potential of N2, O2, Ar + 8%H2 and CO2 nanobubbles in deionized water as a function of pH.

**Table 1.** Solubility of gas, dielectric constant [37] and isoelectric point of nanobubbles determined from the results shown in Figure 1.


The natural pH of CO2 nanobubbles in the deionized water was around 4.2. The IEPs of nanobubbles in the deionized water containing air, Ar + 8%H2, O2, N2, and CO2 had a pH of 4.2, 4.6, 4.7, 5.2 and 5.7, respectively (Table 1). They were determined from the two pH values of zeta potential, assuming a linear relationship. Among all the gases studied, the IEP of CO2 was the highest (i.e., 5.7). This can be explained by the dissolution of the bubble surface, as discussed above, and because the IEP increased with the dielectric constant of gas increase, as listed in Table 1.

#### *3.2. Time Effect of Mean Size and Zeta Potential of Different Gas Nanobubbles and pH and Eh of Nanobubble Suspensions*

The effect of time after the preparation of the nanobubble water on the nanobubble mean diameter of N2, O2, Ar + 8%H2, air and CO2 gas in deionized water (A), 1 mM NaCl aqueous solution (B), 1 mM CaCl2 aqueous solution (C) and 1 mM AlCl3 aqueous solution (D) are shown in Figure 3. The zeta potential of the nanobubbles and the pH and Eh as a function of time are shown in Figures 4 and 5, respectively. In this section, we introduce and discuss the results in the absence of salt followed by the results in the presence of salt. The Eh value change can be correlated to the H2 nanobubble concentration change, O2 nanobubble existence and CO2 nanobubble dissolution in water.

(**A**) Five kinds of gas bubbles' mean diameter in deionized water (**B**) Five kinds of gas bubbles' mean diameter in 1 mM NaCl

aqueous solution

(**C**) Five kinds of gas bubbles' mean diameter in 1 mM CaCl2 aqueous solution

(**D**) Five kinds of gas bubbles' mean diameter in 1 mM AlCl3 aqueous solution

**Figure 3.** Effect of time on nanobubble mean diameter of different kind of gases (N2, O2, Ar + 8%H2 air and CO2) in (**A**) deionized water, (**B**) 1 mM NaCl aqueous solution, (**C**)1 mM CaCl2 aqueous solution and (**D**) 1 mM AlCl3 aqueous solution.

(**A**) Five kinds of gas bubbles' zeta potential in deionized water (**B**) Five kinds of gas bubbles' zeta potential in 1 mM NaCl aqueous solution

**Figure 4.** *Cont.*

(**C**) Five kinds of gas bubbles' zeta potential in 1 mM CaCl2 aqueous solution

(**D**) Five kinds of gas bubbles' zeta potential in 1 mM AlCl3 aqueous solution

**Figure 4.** Effect of time on nanobubble zeta potential of different kind of gases (N2, O2, Ar + 8%H2, air and CO2) in (**A**) deionized water, (**B**) 1 mM NaCl aqueous solution, (**C**) 1 mM CaCl2 aqueous solution and (**D**) 1 mM AlCl3 aqueous solution.

(**B**) Five kinds of gas bubbles' pH and Eh in 1 mM NaCl aqueous solution

(**A**) Four kinds of gas bubbles' pH and Eh in deionized water (**A'**) CO2 gas bubbles' pH and Eh in deionized water

(**C**) Five kinds of gas bubbles' pH and Eh in 1 mM CaCl2 aqueous solution

**Figure 5.** *Cont.*

aqueous solution

**Figure 5.** Effect of time on pH (solid line) and Eh (dashed line) of aqueous solution containing nanobubble of different kind of gases (N2, O2, Ar + 8%H2, air and CO2) in (**A**,**A'**) deionized water, (**B**) 1 mM NaCl aqueous solution, (**C**) 1 mM CaCl2 aqueous solution and (**D**) 1 mM AlCl3 aqueous solution.

#### 3.2.1. Nanobubble Characteristics Change as a Function of Time in Deionized Water

In the absence of salt, the initial mean diameters of the N2, O2, Ar + 8%H2 and air nanobubbles were around 200 nm and gradually increased to 400–530 nm in 30 days. After 30 days, the mean diameter slightly decreased and was stable at 330–480 nm in 60 days (Figure 3A). The zeta potentials of those gases were between −12 and −32 mV for 60 days (Figure 4A). Variation trends of nanobubbles are similar to previous reports. Ulatowski et al. (2019) reported the O2 nanobubbles size produced by porous tube type was 200 nm on the initial day and unchanged after 10 days, while the N2 nanobubble mean size was 300 nm at −11 mV of zeta potential on the initial day and 400 nm at about −12 mV after 35 days [3]. Meegoda et al. (2019) reported the O2 nanobubble mean size produced by the cavitation method was 180 nm at −20 mV of zeta potential on the initial day and the O2 bubble mean size increased to 285 nm at −16 mV after 7 days [9]. Michailidi et al. (2020) reported that with the nanobubbles prepared by the high shear cavitation method, initially the O2 nanobubble mean size was 350 nm and the air nanobubble mean size was 430 nm, while the O2 and air nanobubble sizes increased to 560 and 500 nm, respectively, after two weeks; both bubbles existed at 640 nm size after 2 and 3 months [4]. Although the nanobubble size depends on the preparation methods, the N2 and O2 nanobubbles' size increase with time was similar to our results that the prepared nanobubble sizes of various gases were about 200 nm initially and the nanobubble size increased to 300 to 500 nm with 2 to 4 weeks. Nirmalkar et al. (2018) reported that the number density of air nanobubbles prepared by the acoustic cavitation method gradually decreased; however, about 90 nm of nanobubbles existed for almost one year [6].

On the other hand, the CO2 nanobubble size gradually increased with time and the size changed from 180 nm initially to 350 nm in 5 days (Figure 3A). After 5 days, the CO2 nanobubbles were not stable and could not be detected by DLS because the CO2 dissolved in the water and the electric double layer was compressed; therefore, the bubbles were coagulated due to the Van der Waals interaction (and the hydrophobic interaction) as discussed in the following Section 3.3.1. Ushikubo et al. (2010) reported that with the nanobubbles prepared by the pressurized type, the zeta potential of CO2 nanobubbles was negative at pH 4, initially; however, the bubbles disappeared at pH 6.1 due to the dissociation of CO2 [39]. In our experiment, the zeta potential of CO2 nanobubbles was +9 mV at pH 4.2 on the first day and it gradually increased with time to +17 mV at pH 4.6 after 5 days (Figure 4A). As shown in Figure 5A, the Eh of Ar + 8%H2 nanobubble suspension increased from 250 to 400 mV in one day and also the Eh of N2 nanobubble suspension increased from 330 to 420 mV in one day. The Eh of O2 and air nanobubble suspensions were between 400 and 450 mV initially, then to 60 days. The Eh of the CO2 bubble suspension increased as the pH increased with time (Figure 5A').

#### 3.2.2. Nanobubble Characteristics Change as a Function of Time in Salt Aqueous Solutions

In the presence of salt, the nanobubbles were prepared in 1 mM NaCl, CaCl2 or AlCl3 aqueous solution. In 1 mM NaCl aqueous solution, N2 and air nanobubble sizes were initially 200 nm and gradually increased; they were not observed after 7 days. Meanwhile, O2, Ar + 8%H2, and CO2 bubble sizes were initially 200 to 300 nm, gradually increased and they were not detected after 14 days (Figure 3B). The existence period of CO2 bubbles became longer in the 1 mM NaCl aqueous solution than in the absence of salt (i.e., 5 days in deionized water) (Figure 3A). The zeta potentials of N2, O2, Ar + 8%H2 and air nanobubbles were small, within ±5 mV at initial days to 14 days (Figure 4B), while the zeta potential of CO2 nanobubbles was higher than +12 mV initially and was kept +10 mV for 14 days (Figure 4B). Ke et al. (2019) reported that with the nanobubbles produced by the pressurized type, the initial N2 nanobubble average size was about 200 nm in the 0.1 mM NaCl aqueous solution [2] and the size was similar to our result in the 1 mM NaCl condition (Figure 3B). In 1 mM NaCl aqueous solution, Leroy et al. (2012) measured that the zeta potential of H2 gas was −20 mV at pH 6 [40] and Yang et al. (2001) reported the initial zeta potential of H2 nanobubbles was −30 mV [41]. In our measurement of Ar + 8%H2 nanobubble, the zeta potential was much smaller at −5 mV at pH 6 in 1 mM NaCl (Figure 3B). Meegoda et al. (2019) reported that the O2 nanobubble mean size prepared by the high shear cavitation type was 214 nm at −22 mV of zeta potential initially and the O2 bubble mean size was almost the same at 219 nm at −14 mV in 7 days [9]. In our experiment, the O2 bubble mean size was about 300 nm for 14 days (Figure 3B), similar to the above-mentioned previous result; however, the zeta potential was +6 mV (Figure 4B). The Eh value increased from 220 to 360 mV with the Ar + 8%H2 nanobubble suspension in one day, which was the same as nanobubbles in the deionized water suspension. The Eh of the N2 nanobubble suspension was 360 mV and the Eh of other nanobubble suspensions was a little high at about 450 mV in 7 days (Figure 5B). The pH of CO2 was from 4.3 to 4.4 and the Eh, pH, zeta potential and mean size of the CO2 suspension were almost constant in the 1 mM NaCl aqueous solution compared with its suspension in deionized water.

In 1 mM CaCl2 salt aqueous solution, O2 was not stable after 7 days, while N2, Ar + 8%H2 and air bubbles could be observed until 14 days, but no more were detected after 14 days (Figure 3C). The CO2 nanobubbles existed until about 20 days in 1 mM CaCl2. The absolute value of the zeta potential of N2, O2, Ar + 8%H2 and air nanobubbles was low, at less than 10 mV, while the zeta potential of CO2 was higher than +10 mV (Figure 4C). Cho et al. (2005) reported that the air nanobubble size prepared by the sonicating method was 850 nm at −8 mV zeta potential [42], and their reported results were similar to our results after one day (Figures 3C and 4C). The size of the air, N2 and O2 nanobubbles in 1 mM CaCl2 salt aqueous solution increased extremely with time. Yang et al. (2001) reported that the zeta potential of the H2 nanobubble was −5.5 mV [8], and their result was similar to our result of Ar + 8%H2 nanobubble after one day (Figure 4C). For the Eh value (Figure 5C), the Eh of the 8%H2+Ar nanobubble suspension increased from 190 to 370 mV in one day, same as the nanobubbles in the other water suspensions. The Eh of N2, O2, air and CO2 nanobubble suspensions were between 400 and 450 mV and the pH of N2, O2, Ar + 8%H2 and air nanobubble suspensions were between 5.5 and 6. The pH of CO2 solution increased from 4.4, initially, to 5.3 after 7 days (Figure 5C).

In 1 mM AlCl3 salt aqueous solution, N2 and O2 were not stable after 7 days, while air, CO2 and Ar + 8%H2 nanobubbles could be observed until 14 days; however, they were not detected after 14 days. The size of the Ar + 8%H2 nanobubble was more stable, at between 200 and 300 nm for 14 days, compared with other bubbles (Figure 3D). The zeta potentials of all prepared bubbles were positive (Figure 4D) possibly due to the presence of high valence ions (Al3+) adsorbing on the bubble surfaces. Yang et al. (2001) reported that

the zeta potential of H2 nanobubble was +12 mV [41], and their result was similar to our result of +15 mV of Ar + 8%H2 nanobubbles at the initial day (Figure 4D). Han et al. (2006) reported that the zeta potential of O2 nanobubbles prepared by electrolysis was +20 mV at pH 6 in 10 mM NaCl + 1 mM AlCl3 aqueous solution [43]. The trivalent ion Al3+ ion caused the positive charge in all our prepared nanobubbles (Figure 4D). For the Eh value (Figure 5D), the Eh of the Ar + 8%H2 nanobubble suspension increased from 190 to 410 mV in one day, same as the gas nanobubbles in the other water suspensions. The Eh of N2, O2, air and CO2 nanobubble suspensions were between 450 and 500 mV, and were higher than the Eh of the Ar + 8%H2 suspension. The pH of N2, O2, Ar + 8%H2 and air were from 4.15 to 4.35, which were lower than the pH of other suspensions. The pH of CO2 slightly increased from 4.0 to 4.3 in 5 days (Figure 5A') and its increase was limited compared with other suspensions in the salt aqueous solution and deionized water (Figure 5B).

The stable days of N2, O2, 8% H2+Ar, air and CO2 nanobubbles produced in different solutions are listed in Table 2. The N2, O2, Ar + 8%H2 and air nanobubbles in 1 mM of the three salts with different valences (+1, +2 and +3) decreased the existence period for 1 to 2 weeks compared with deionized water in the absence of salt (>60 days). On the other hand, CO2 bubbles existed only 5 days in deionized water in the absence of salt; however, CO2 bubbles existed 14 days in the 1 mM salt aqueous solution. The stability differences will be further discussed in the following section.

**Table 2.** Stable days of gas nanobubbles produced in different solutions. The minimum days are given.


#### *3.3. Nanobubble Stabilization Estimation by Using the Extended DLVO Theory*

The thickness of the electric double layer (Debye length = 1/*κ*) of the nanobubble is calculated in the next formula:

$$\kappa = \sqrt{\frac{2nz^2e^2}{\varepsilon\_r \varepsilon\_0 kT}}\tag{9}$$

where *κ* is the Debye–Hückel parameter, *n* is the concentration of anion or cation in the solution and is equal to 1000 *NAC* (*NA* is the Avogadro number and *C* is concentration mol/L), *z* is the valence of ion, *e* is the electron charge, *ε<sup>r</sup>* is the relative dielectric constant and *ε*<sup>0</sup> is the permittivity of vacuum.

The calculated Debye lengths of the nanobubbles studied in our experiment are listed in Table 3 in order to discuss differences in the nanobubble stability and to have some idea of how to keep the nanobubble stable for a longer duration. The Debye length is 300 nm at pH 6 around N2, O2 and Ar + 8%H2 nanobubbles in deionized water, and 140 nm at pH 4.2 and 70 nm at pH 4.7 around CO2 nanobubbles. Since the characteristics of CO2 nanobubbles are different from others (Figures 3–5), the Debye length (1/*κ*) around CO2 bubbles and HCO3 −, OH−, CO3 <sup>2</sup><sup>−</sup> ion concentration in the presence of CO2 nanobubbles in deionized water as a function of pH were also calculated and shown in Figure 6 in order to discuss their relationship. The concentration of various ions was calculated from using Equations (3)–(8).


**Table 3.** Debye length (1/*κ*, nm) of nanobubbles prepared in aqueous solutions at natural pH.

**Figure 6.** Debye length (1/*κ*) around CO2 bubbles (left axis) and HCO3 −, OH−, CO3 <sup>2</sup><sup>−</sup> ion concentration in the presence of CO2 nanobubbles in deionized water (right axis) as a function of pH.

A thick electric double layer (300 nm, Table 3) is present around the N2, O2 and Ar + 8%H2 nanobubbles in deionized water and can stabilize those nanobubbles with high enough zeta potential (i.e., −16 to −32 mV, Figure 4A). On the other hand, the Debye length of CO2 nanobubbles (70 nm at pH 4.7 in DI water, Table 3) and other nanobubbles in salt aqueous solution (3–10 nm, Table 3) are thin and thus, the influence of the Van der Waals and hydrophobic attraction can be more dominant in coagulating and coalescing the nanobubbles. In order to quantify the influence of the electrical double layer potential, the Van der Waals potential and the hydrophobic interaction potential, their potential energies were calculated by using the extended DLVO theory as described in the following sections. Tchaliovska et al. investigated the thickness of thin flat foam films formed in aqueous dodecyl ammonium chloride solution [44]. Angarska et al. showed that the value of the critical thickness of the foam film was sensitive to the magnitude of the attractive surface forces acting on the film using sodium dodecyl sulfate and the total disjoining pressure operative on the films, which was expressed as a sum of the Van der Waals and hydrophobic contributions [45].

Figure 7 shows the two nanobubble positions and geometries considered in our potential calculation. The total potential energy *VT* between two nanobubbles can be given by the potential energy due to the Van der Waals interaction (*VA*), hydrophobic interaction (*Vh*) and the electrostatic interaction (*VR*) as follows [46–50];

$$V\_T = V\_A + V\_h + V\_{R\prime} \frac{V\_T}{kT} = \frac{V\_A + V\_h + V\_R}{kT} \tag{10}$$

*VA* + *Vh* is shown in the follow equation:

$$V\_A + V\_h = -\frac{A+K}{6} \left[ \frac{2a\_1a\_2}{R^2 - (a\_1+a\_2)^2} + \frac{2a\_1a\_2}{R^2 - (a\_1-a\_2)^2} + \ln \left( \frac{R^2 - (a\_1+a\_2)^2}{R^2 - (a\_1-a\_2)^2} \right) \right] \tag{11}$$

where the Hamaker constant *<sup>A</sup>* for air in water whose value of air–water–air of 3.7 × <sup>10</sup>−<sup>20</sup> J [46] was used to investigate our system in nanobubble–water–nanobubble. The *K* is a hydrophobic constant for air in water. Yoon and Aksoy calculated the *K* using the extended DLVO theory [49] while Wang and Yoon measured the *K* as a function of the SDS concentration at different kinds of NaCl concentrations and showed that *K* was estimated at 10−<sup>17</sup> J in the absence of salt and 10−<sup>19</sup> J in the 1 mM NaCl aqueous solution [50]. In our calculation, 10−<sup>19</sup> J of *K* was used in 1 mM CaCl2 and AlCl3 aqueous solutions and 10−<sup>18</sup> J of *K* was used in deionized water containing CO2 nanobubbles by considering the reference values [50].

**Figure 7.** Position and geometry of two nanobubbles.

The radius of nanobubbles is *a*<sup>1</sup> and *a*<sup>2</sup> and their distance *R* = *a*<sup>1</sup> + *a*<sup>2</sup> + *h* is defined as shown in Figure 7, and *h* is the surface-to-surface distance between two nanobubbles. The potential energy of *VR* is expressed as follows,

$$V\_R = -\frac{\pi \varepsilon\_r \varepsilon\_0 a\_1 a\_2 \left(\psi\_1^2 + \psi\_2^2\right)}{a\_1 + a\_2} \left[\frac{2\psi\_1 \psi\_2}{\psi\_1^2 + \psi\_2^2} \ln \frac{1 + \exp(-\kappa h)}{1 - \exp(-\kappa h)} + \ln\{1 - \exp(-2\kappa h)\}\right] \tag{12}$$

where *ψ*<sup>1</sup> *and ψ*<sup>2</sup> are surface potentials of the nanobubbles of radii *a*<sup>1</sup> and *a*2, respectively.

#### 3.3.1. Stability Calculation of Nanobubbles in Deionized Water

The mean diameters of the prepared nanobubbles of N2, O2, Ar + 8%H2 and air in deionized water were between 170 and 230 nm, as shown in Figure 3. With time, these diameters gradually increased, and they were between 390 and 530 nm after 30 days. Then, they slightly decreased to between 330 and 480 nm and were stable after 60 days. As the mean bubble diameter was measured in our experiment, the increase in bubble size indicates that the small size nanobubbles coalesce with other bubbles, and those small bubbles disappear.

The above-mentioned phenomenon is well known as the Ostwald ripening effect [51], which explains the deposition of a smaller object on a larger object due to the latter being more energetically stable than the former. By using Lemlich's theory [52], it can be explained by using next equation [53]:

$$\frac{da}{dt} = \mathcal{K}' \left( \frac{1}{p} - \frac{1}{a} \right) \tag{13}$$

where *a* is the bubble radius, *t* is time, *K* is a constant and *p* is instantaneous bubble mean radius. This equation tells us that the bubbles with a radius larger than *a* grow, and those with a radius smaller than *p* shrink. By using the results shown in Figure 2 (zeta potential) and Figure 3 (bubble size), the bubble stability is further discussed by calculating the total potential energy between two nanobubbles.

Figure 8 shows the potential energies between two nanobubbles in deionized water. Figure 8A shows the potential energies between two same-size gas bubbles of 200 nm of

Ar + 8%H2 on the initial day. The total potential energy *VT* is 20 kT at 300 nm (the Debye length) and the maximum potential is 30 kT. The other bubbles of air, N2 and O2 bubble show almost same potential curves. It indicates that the two bubbles have enough of a potential barrier to disperse each other. Figure 8B shows the potential energies between 450 nm of two same-size bubbles of Ar + 8%H2 after 30 days. The maximum total potential energy *VT* is 15 kT, i.e., the threshold total potential energy to determine coagulation or dispersion [54], and it indicates that Ar + 8%H2 nanobubbles would be stable, although the absolute value of zeta potential is smallest in Ar + 8%H2 nanobubbles (−13 mV, Figure 4A) compared with the ones of other O2, N2 and air nanobubbles bubbles, shown in Figure 4.

(**A**) Two 200 nm Ar + 8%H2 gas bubbles on the prepared day (**B**) Two 450 nm Ar + 8%H2 gas bubbles after 30 days

**Figure 8.** Total potential energy *VT* as a function of the distance between two nanobubbles at −24 mV zeta potential on the prepared day for Ar + 8%H2 (**A**) and −12 mV after 30 days for Ar + 8%H2 nanobubbles (**B**) in deionized water. (**A**) Two 200 Ar + 8%H2 gas bubbles on the prepared day, (**B**) two 450 nm Ar + 8%H2 gas bubbles after 30 days.

> When the CO2 nanobubbles were prepared, their mean size was 160 nm and zeta potential was low, i.e., +9 mV. Figure 9A shows the potential energies between two 160 nm bubbles. The total potential barrier was not appeared at 140 nm the Debye length. Figure 9B shows the potential energies between two 350 nm bubbles after 5 days. As the total potential barrier was also not appeared, these larger bubbles would also coagulate, coalesce and increase the size and, thus, CO2 bubbles could not be observed after 5 days.

(**A**) Two 160 nm CO2 gas bubbles on the prepared day (**B**) Two 350 nm CO2 gas bubbles after 5 days

**Figure 9.** Total potential energy *VT* as a function of the distance between two nanobubbles at +9 mV zeta potential on the prepared day for CO2 (**A**) and +5 mV after 5 days for CO2 nanobubble (**B**) in deionized water. (**A**) Two 160 nm bubbles on the prepared day, (**B**) two 350 nm bubbles after 5 days.

#### 3.3.2. Stability Calculation of Nanobubbles in Salt Aqueous Solutions

In the 1 mM NaCl aqueous solution, the absolute values of the zeta potential of O2, N2, Ar + 8%H2 and air nanobubbles were less than 10 mV; the total potential energy barrier could not appear. Figure 10A shows the potential energies between two 160 nm N2 nanobubbles. The total potential energy *VT* is negative at any distance and does not show the potential barrier; therefore, the nanobubbles are not stable and can coagulate. Figure 10B shows the potential energies between two 200 nm CO2 nanobubbles. Although the absolute zeta potential was slightly larger than the sN2 nanobubbles at the Debye length 10 nm, there is no potential barrier indicating the attraction between the two bubbles. The CO2 nanobubble size existed between 200 and 400 nm for 14 days, as shown in Figure 3B; however, after 14 days, the CO2 bubble was disappeared.

**Figure 10.** Total potential energy *VT* as a function of the distance between two nanobubbles at −5 mV zeta potential on the prepared day for N2 (**A**) and +12 mV on the prepared day for CO2 nanobubble (**B**) in 1 mM NaCl aqueous solution. (**A**) Two 160 nm N2 bubbles on the prepared day, (**B**) two 200 nm CO2 bubbles on the prepared day.

> In the 1 mM CaCl2 aqueous solution, as the absolute values of the zeta potential of O2, N2, Ar + 8%H2 and air nanobubbles were less than 10 mV (Figure 4C), like the one in 1 the mM NaCl aqueous solution (Figure 4B), the total potential energy barrier did not appear. Figure 11A shows the potential energies between 230 nm of two same-size N2 nanobubbles on the prepared day. The total potential energy *VT* was negative at any distance and did not show the potential barrier, and thus it indicates that the nanobubbles were not stable due to their coagulation. The zeta potential of CO2 was more than +10 mV, as shown in Figure 4C; however, the potential barrier still did not appear because the Debye length around the CO2 bubbles decreased in the presence of 1 mM CaCl2 salt (i.e., 5 nm vs. 40 nm (pH 4.2) in deionized water, Table 2, while the CO2 nanobubble mean size 1 mM CaCl2 was stable at about 300 nm until 14 days (Figure 3C), same as the one in the 1 mM of NaCl aqueous solution (Figure 3B). Figure 11C shows the potential energies between two 750 nm N2 nanobubbles after 3 days in a 1 mM CaCl2 aqueous solution. With the smaller zeta potential −7 mV, the potential barrier is not appeared, and it corresponds to the further coagulation/coalescence of the bubbles (Figure 3C).

> In the 1 mM AlCl3 aqueous solution, Figure 11B shows the potential energy between two 170 nm Ar + 8%H2 nanobubbles on the prepared day. Although the absolute value of the zeta potential was higher (+16 mV, Figure 4D) than the one of Ar + 8%H2 nanobubbles prepared in a CaCl2 aqueous solution (−2 mV, Figure 4C), there was no potential barrier due to the small electrostatic interaction potential (<15 kT) and the bubbles may be unstable. In Figure 11D, with the N2 nanobubble size of 300 nm and the higher zeta potential (+25 mV, Figure 4D), the potential barrier was not also appeared.

(**A**) Two 230 nm N2 gas bubbles at ƺ10 mV in 1 mM CaCl2 aqueous solution on the prepared day (**B**) Two 170 nm Ar + 8%H2 gas bubbles at +16 mV in AlCl3 aqueous solution on the prepared day

(**C**) Two 750 nm N2 gas bubbles at -7mV in 1 mM CaCl2 aqueous solution at 3rd day

**Figure 11.** Total potential energy *VT* depending on the distance between two nanobubbles in the salt aqueous solution. (**A**) Two 230. nm N2 nanobubbles at −10mV in 1 mM CaCl2 aqueous solution on the prepared day, (**B**) Two 170 nm Ar + 8%H2 nanobubbles at +16 mV in AlCl3 aqueous solution on the prepared day, (**C**) Two 750 nm N2 gas bubbles at −7mV in 1 mM CaCl2 aqueous solution at 3rd day, (**D**) Two 300 nm N2 nanobubbles at +25 mV in 1 mM AlCl3 aqueous solution after 3 days.

> In the case of N2, O2, Ar + 8%H2 and air in 1 mM salt concentration, the presence of Ca2+ ion from CaCl2 salt decreased the absolute value of the zeta potential by its adsorption on the bubble surface and a stronger coagulation could happen by thinning the electric double layer at the natural pH (without pH adjustment). Here, the overall bubble stability phenomena are summarized. Except for the 1 mM Ca2+ ion aqueous solution, the stability phenomena nearly followed the Shultz–Hardy rule. Comparing bubble size and existence time, the stability order of the N2, O2, Ar + 8%H2 and air nanobubbles identified from our study was: no salt addition > 1 mM NaCl > 1 mM AlCl3 > 1 mM CaCl2 in deionized water. On the other hand, the stability order of CO2 nanobubbles was: 1 mM NaCl > 1 mM CaCl2 > 1 mM AlCl3 > no salt addition, due to the effect of CO2 dissolution in deionized water, as discussed in Section 3.3.

#### 3.3.3. Nanobubble Movement in Salt Aqueous Solution

The nanobubbles move by Brownian motion. The Brownian diffusion can be described in the following Equations (14) and (15) where the bubble movement distance is Δ*x*<sup>2</sup> and Equation (1) is incorporated to define *DT*.

$$
\overline{\Delta \mathbf{x}^2} = 2D\_T \mathbf{t} = 2\mathbf{t} \mathbf{k} \mathbf{T} / 3\pi \mathbf{\tau} \eta \mathbf{d}\_h \tag{14}
$$

$$
\sqrt{\overline{\Delta x^2}} = \sqrt{\frac{2tkT}{3\pi\eta d\_h}}\tag{15}
$$

If the bubble size *dh* is small, the movement distance becomes longer since they have inverse correlation, as described in Equations (14) and (15).

On the other hand, a large, coalesced bubble floats, following the terminal velocity *u* defined by the Stokes equation of laminar flow:

$$
\mu = \frac{d\_h^2 \rho g}{18\eta} \tag{16}
$$

where *ρ* is the density of water and *g* is the gravitational acceleration.

Bubbles experiencing diffusion change their movement direction frequently in a short time. On the other hand, the bubble displacement due to the buoyancy force is always pointing upwards against the gravitational force. The terminal velocity depends on the bubble diameter, as shown in Figure 12, based on our calculation using Equation (16). A bubble smaller than 1 <sup>μ</sup>m can experience very low terminal velocity (<1 × <sup>10</sup>−<sup>6</sup> m/s) which prevents bubbles from floating against the buoyancy force. Our results partially agree with the previous literature. Nirmalkar et al. (2018) reported that the number density of about a 100 nm mean size of nanobubbles gradually decreased over one year, and still existed at about 100 nm after one year [8].

**Figure 12.** Relationship between the terminal velocity of a floating bubble and bubble size.

Figure 13 shows the schematic diagrams of nanobubble stability in deionized water (Figure 13A) and the nanobubble size change in a salt aqueous solution (Figure 13B). The nanobubbles in deionized water stably exist for two months by the Brownian motion. On the other hand, in the nanobubbles in a 1 mM salt aqueous solution, the nanobubbles coalesced with each other and increased the size, and the large-size bubbles gradually floated and disappeared after one or two weeks.

The surface charge on the front direction of a floating larger bubble can decrease, while the tail direction of the bubble retains more ions; thus, the electric double layer charge distribution can be distorted [55]. The bubbles in the direction of a floating large bubble can interact with that large bubble. If there is a small nanobubble in the front direction of a moving bubble, it can coagulate and coalesce with large bubbles and the coalesced large bubbles can be disappeared. This phenomenon occurred for the larger size of the bubbles in a salt aqueous solution with time, as observed in this study.

**Figure 13.** Schematic diagrams of nanobubble stability in deionized water (**A**) and nanobubble size change in salt aqueous solution (**B**) depending on time.

#### **4. Conclusions**

This research investigated five kinds of gas nanobubbles (N2, O2, Ar + 8%H2, air and CO2) for their long-term stability, to consider the application of the nanobubbles containing aqueous solutions. Each gas in a cylinder was injected into deionized water or a 1 mM of salt aqueous solution, and nanobubbles were prepared by the hydrodynamic cavitation method. The mean diameter, zeta potential of bubbles, pH and Eh of nanobubble suspensions were measured and these characteristics changed, as the function of time was also studied.


**Author Contributions:** Y.Z.: Conceptualization, data curation, writing; Z.H.: Data curation; C.H.: Funding acquisition, Investigation; K.W.: Formal analysis; Y.W. (Youbin Wang): Formal analysis; N.L. and Q.F.: Methodology; G.D. and A.O.: Writing and editing; Y.W. (Yuezhou Wei): Supervision, investigation; T.F.: Conceptualization, supervision, writing, reviewing and editing. All authors have read and agreed to the published version of the manuscript.

**Funding:** One part of this research was funded by the Natural Science Foundation of China, grant number NSFC 21976039.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The authors confirm that the data supporting the findings of this study are available within the article.

**Acknowledgments:** The authors appreciate the reviewers.

**Conflicts of Interest:** The authors declare no competing interests.

#### **Nomenclature**



#### **References**


### *Article* **The Influence of Air Nanobubbles on Controlling the Synthesis of Calcium Carbonate Crystals**

**Yongxiang Wu 1, Minyi Huang 1, Chunlin He 1, Kaituo Wang 1, Nguyen Thi Hong Nhung 1, Siming Lu 2, Gjergj Dodbiba 3, Akira Otsuki 4,5 and Toyohisa Fujita 1,\***


**Abstract:** Numerous approaches have been developed to control the crystalline and morphology of calcium carbonate. In this paper, nanobubbles were studied as a novel aid for the structure transition from vaterite to calcite. The vaterite particles turned into calcite (100%) in deionized water containing nanobubbles generated by high-speed shearing after 4 h, in comparison to a mixture of vaterite (33.6%) and calcite (66.3%) by the reaction in the deionized water in the absence of nanobubbles. The nanobubbles can coagulate with calcite based on the potential energy calculated and confirmed by the extended DLVO (Derjaguin–Landau–Verwey–Overbeek) theory. According to the nanobubble bridging capillary force, nanobubbles were identified as the binder in strengthening the coagulation between calcite and vaterite and accelerated the transformation from vaterite to calcite.

**Keywords:** nanobubbles; calcium carbonate; transformation; crystal-control; extended DLVO theory; capillary force

#### **1. Introduction**

Nanobubbles are described as a gaseous domain in the liquid phase with a diameter of less than 1 μm [1] with unusual longevity lasting for several weeks or even months, which contradicts the classical Epstein–Plesset theory prediction that nanobubbles cannot exist and should disappear in a few milliseconds or microseconds [2,3]. Until now, significant progress has been made on the stability of nanobubbles [2–4], which has attracted the attention of experimentalists and theorists. There are more and more studies involving possible explanations for the stability of nanobubbles that have been proposed [5]. Since it was discovered that the surface of nanobubbles had a high magnitude of zeta potential, the charge stabilization model has aroused the attention of researchers [6–8]. In recent years, the mechanism of charge stabilization appears to be the most reasonable to rationalize nanobubble stabilization. An electrical double layer is formed on the surface of the nanobubble in water containing free moving ions, which serves as the main energy supply to withstand the coalescence of multiple bubbles and generates a repulsive pressure to balance the Laplace pressure [2].

The majority of nanobubble research has recently been rapidly increasing in many kinds of scientific disciplines. The unique physicochemical properties of the nanobubbles, such as long-term stability [2,3,6], high zeta potential [3,7], generation, and degradation of free radicals when collapsing [9], has aroused interest in many areas. Therefore, numerous reports have been published on the application of nanobubbles, including medicine [10,11], environment [9,12], flotation [13,14], and materials [15]. In medicine, nanobubbles have

**Citation:** Wu, Y.; Huang, M.; He, C.; Wang, K.; Nhung, N.T.H.; Lu, S.; Dodbiba, G.; Otsuki, A.; Fujita, T. The Influence of Air Nanobubbles on Controlling the Synthesis of Calcium Carbonate Crystals. *Materials* **2022**, *15*, 7437. https://doi.org/10.3390/ ma15217437

Academic Editor: Carlos Javier Duran-Valle

Received: 2 October 2022 Accepted: 20 October 2022 Published: 23 October 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/).

been used to protect proteins from surface-mediated denaturation [16] as well as for ultrasound imaging and intracellular drug delivery [17]. H2 nanobubbles in water were discovered to inhibit tumor cell growth [18]. The application of nanobubbles in environmental studies is well-known, such as the removal of heavy metal ions [19], wastewater treatment [20], removal of pollutants from water [21]. As mentioned above, nanobubbles are thought to have a wide range of applications. However, the application of nanobubbles in material science and engineering is uncommon, necessitating further exploration of their potential applications.

To some extent, nanobubbles can change the state and properties of water, as mentioned in many papers. Hydrogen nanobubbles in aqueous solution have a very low Eh (oxidation and reduction potential), which can be used to prevent oxidation [8]. Because a large number of hydroxyl radicals are generated when oxygen nanobubbles in aqueous solution break, they have a high oxidation state [8]. Due to their unique physicochemical properties, in this research, nanobubbles in water were used as a reaction solvent instead of pure water during calcium carbonate (CaCO3) synthesis by the double decomposition method to study the influence of nanobubbles on the growth of crystalline CaCO3.

Calcium carbonate (CaCO3), one of nature's most common and widely dispersed materials, is an important component of the global carbon cycle [22]. Calcium carbonate has three crystalline phases vaterite, aragonite, and calcite [23]. In these crystals, vaterite is the least stable phase in the water, and can spontaneously change into calcite [24]. Calcite is the most common form because it is thermodynamically stable [23]. The formation of vaterite to calcite can be influenced by reaction parameters, such as supersaturation, temperature, reaction time, and additive use [25].

Furthermore, calcium carbonate has a wide range of applications that depend on its shape and morphology. Therefore, it is important to control the crystal shape or morphology of CaCO3 [26]. Aragonite whiskers with high aspect ratios, for example, have been in great demand for the reinforcement of polymer materials [26]. Flexible and deformable calcium carbonate is also extensively used as a filler and pigment in papermaking [27]. Surface-functionalized calcium carbonate particles can also be adapted to create novel catalytic materials [28]. Non-spherical vaterite particles are appealing as solid supports for targeted and extended drug delivery [29]. In wastewater treatment, vaterite has been explored as it is excellent at removing several heavy metal ions [30–32]. Porous vaterite and cubic calcite aggregated calcium carbonate are used for Cu2+ heavy metal removal [32]. Vaterite particles for strontium removal were also investigated and reported [31]. However, Sasamoto reported that porous calcite has potential as an adsorbent with a fast reaction rate and large adsorption amount in the removal of heavy metal ions [33].

There are several preparation methods for CaCO3, including the carbonation method [34], double decomposition method [35], and the thermal decomposition of calcium bicarbonate [36]. Among them, the double decomposition method is a common way to produce CaCO3. In this method, calcium chloride (CaCl2) is exploited as a calcium source, and sodium carbonate (Na2CO3) is used as a carbon source [35]. In some papers, different kinds of water-soluble additives have been used to control the crystal formation of the CaCO3, such as sucrose [37], anion surfactants [38], sodium dodecylbenzene sulfonate [39] and para-aminobenzoic acid [40]. Nanobubbles are clean and environmentally friendly. On the other hand, there is little research on the effect of nanobubbles on the crystal formation of calcium carbonate.

Due to the stability of the high zeta potential and the nanobubble bridging capillary force, nanobubbles could affect the crystal transition of calcium carbonate. According to previous studies [8], nanobubbles containing different gases, except CO2, had similar physical properties. Therefore, only air nanobubbles were applied in the present study. In future studies, CO2 and other gases would be discussed.

In this paper, air nanobubbles were produced using a hydraulic cavitation method [41,42]. The double decomposition method was used to produce calcium carbonate. To control the formation of CaCO3 crystals, nanobubbles were used instead of traditional additives

in this research. The DLVO (Derjaguin–Landau–Verwey–Overbeek) theory [43,44] was used to evaluate the stability of the colloidal CaCO3 particles containing electrolytes in aqueous solution, while extended DLVO theory incorporating hydrophobic interaction energy [8] was used to evaluate the stability of the nanobubbles and the interaction between nanobubbles and CaCO3, because of the hydrophobic surface property of nanobubbles. According to the experimental results and calculations, the possible kinetics and mechanisms of nanobubbles' influence on the transformation were proposed in detail.

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

#### *2.1. Materials*

Deionized water with a resistivity of 18.2 MΩ·cm−<sup>1</sup> prepared by the Classic Water Purification System from Hitech instruments Co., Ltd. (Shanghai, China) was used for all the experiments. Calcium chloride (CaCl2, AR, 96.6%, Guangdong Guanghua Sci-Tech Co., Ltd., Shantou, China), and anhydrous sodium carbonate (Na2CO3, AR, 99.8%, Guangdong Guanghua Sci-Tech Co., Ltd., Shantou, China) were obtained from Guangdong Guanghua Sci-Tech Co., Ltd. (Shantou, China) and were of analytical purity. The pH of the aqueous solution was modified by the addition of hydrochloric acid (HCl, AR, 36%~38%, CHRON CHEMICALS, Chengdu, China) or sodium hydroxide (NaOH, AR, 96%, Guangdong Guanghua Sci-Tech Co., Ltd., Shantou, China).

To investigate the zeta potential of the vaterite and calcite, 100% of the calcite prepared in this work was used. 94% of the artificial vaterite was prepared by bubbling CO2 into 0.05 M CaCl2 aqueous solution, and ammonia (NH4OH, 30%) was used to adjust the solution pH to 10 while the solution mixture was vigorously stirred in a beaker continuously. After 15 min, the particles were rinsed with deionized water and filtered twice, and then dried in an oven at 110 ◦C for 0.5 h and then at 70 ◦C for 12 h.

The reaction of vaterite synthesis is as follows [34]:

$$\text{CaCl}\_2 + \text{CO}\_2 + 2\text{ NH}\_4\text{OH} \rightarrow \text{CaCO}\_3 + 2\text{ NH}\_4\text{Cl} + \text{H}\_2\text{O} \text{ (pH} = 10\text{)}$$

According to the SEM image and the XRD pattern, shown in Figure S1, spherical calcium carbonate was vaterite with high purity (94%).

#### *2.2. Methods*

#### 2.2.1. Preparation of Nanobubble Suspension

The bulk of nanobubbles were generated by high-speed shearing using a purpose-built device (self-made equipment) which consisted of a plastic stent, an electric high-speed motor in a plastic bracket, a glass beaker, and a blade made of stainless steel as shown in Figure S2. The working principle was as follows: high-speed rotation of the blade above the motor drive creates high- and low-pressure areas and forms pressure gradients, which forms cavities in the low-pressure zones where the nanobubbles are produced [41,45].

In a typical process, the device was operated three times by adding 700 mL deionized water at 1 min intervals before producing nanobubbles. The motor frequency was initially set to 20,000 rpm. After this, 700 mL deionized water was poured into the glass beaker. The lid was put on the beaker to prevent the gas from escaping and avoid impurities in the water. The device was operated ten times, with 1 min intervals. Nanobubbles were identified after completing these steps. After generation, the solution was transferred to another glass container and sealed for the next experiment.

#### 2.2.2. Calcium Carbonate Synthesis

The flow chart of calcium carbonate synthesis is shown in Figure 1. Firstly, nanobubbles were produced, CaCl2 and Na2CO3 were then dissolved into the deionized water containing nanobubbles at 20 ◦C. The initial concentrations of CaCl2 and Na2CO3 in deionized water were 0.05 M (M = mol/L). 200 mL Na2CO3 solution (0.05 M) was added to 200 mL CaCl2 solution (0.05 M), and the solution mixture was vigorously stirred in the

beaker continuously. Upon the completion of the designated reaction time (i.e., 0.5, 1, 2, 3, 4 h), the particles were rinsed with deionized water and filtered twice. After the last washing and suction filtration, the sedimented particles were dried in an oven at 110 ◦C for 0.5 h and then at 70 ◦C for 12 h before particle characterization.

**Figure 1.** The flow chart of the synthesis of CaCO3.

#### 2.2.3. Characterization

The size of the nanobubbles was measured by the dynamic light scattering method (DLS, NanoBrook Omni, Brookhaven Instruments, Holtsville, NY, USA) and the concentration of nanobubbles was measured by the particle trajectory method (Nano-Sight, NS300, Malvern). The zeta potential measurements were made with the use of the phase analysis light scattering method (NanoBrook Omni, Brookhaven Instruments). After the synthesized CaCO3 particles were deposited onto a stub and coated with gold, the polymorphs and morphologies were characterized by scanning electron microscopy (SEM, SEM-EDS, Phenom ProX, Netherlands) at an accelerating voltage of 15 kV. The samples obtained in Section 2.2.2 were directly used for XRD measurements. The X-ray diffraction (XRD, Rigaku D/MAX 2500V, Japan's neo-Confucianism) was used to verify the existence of vaterite and calcite and calculate the mass fraction of calcite within the synthesized CaCO3 at Cu Kα target (λ = 0.1542 nm), at 40 kV and 30 mA with a scanning speed of 8◦ 2θ/min. The CaCO3 particles with potassium bromide tableting were measured by Fourier transform infrared spectroscopy (FT-IR, Nicolet iS 50, Thermo Scientific) to verify the existence of vaterite and calcite using the transmission mode with a scanning speed of 20 spectra/s. Typical spectral resolution was 0.25 cm<sup>−</sup>1.

To quantify different crystalline calcium carbonates and calculate the mass fraction of calcite within the synthesized CaCO3, the intensities of the (110), (112), and (114) crystallographic planes for vaterite (I110V), (I112V), and (I114V), and the (104) crystallographic plane for calcite (I104C), the semi-quantitative phase compositions were calculated by using Equations (1) and (2) [46]:

$$\mathbf{X}\_{\rm V} = (\mathbf{I}\_{110\rm v} + \mathbf{I}\_{112\rm v} + \mathbf{I}\_{114\rm v}) / (\mathbf{I}\_{104\rm c} + \mathbf{I}\_{110\rm v} + \mathbf{I}\_{110\rm v} + \mathbf{I}\_{110\rm v}) \tag{1}$$

$$
\chi\_{\mathbf{c}} = 1 - \chi\_{\mathbf{v}} \tag{2}
$$

where Xv is the mass fraction of vaterite, Xc is the mass fraction of calcite, I104C is the intensities of the (104) crystallographic plane for calcite, and I110v, I112v and I114v are the intensities of the (110), (112), and (114) crystallographic planes for vaterite, respectively.

#### **3. Results**

#### *3.1. Generation of Air Nanobubbles*

Because the whole reactor was sealed, impurities were unlikely to be found in this solution. Furthermore, the content of impurities in water was extremely low, indicating

that the colloid particles observed in the solution were nanobubbles. Due to the Tyndall effect based on the light scattering of colloidal sized objects being illuminated with a laser beam after high-speed shearing, a clear light path was observed in the solution, demonstrating that plenty of colloid particles were produced. In this experiment, the hydrodynamic diameter of the air nanobubbles was measured by DLS. The mean diameter of the nanobubbles was about 83.6 nm in the absence of any salt addition. As shown in Table 1, in the 0.05 M CaCl2 and 0.05 M Na2CO3 nanobubble aqueous solutions, the mean diameter was 126.2 nm and 101.8 nm, respectively. The pH of the Na2CO3 aqueous solution was 11.7, while the pH of the deionized water and CaCl2 aqueous solution was about 6. As Ca2+ tended to be adsorb onto the surface of nanobubbles in a CaCl2 aqueous solution, the nanobubbles' zeta potential was positive (7.3 mV). As shown in Figure 2, the concentration of nanobubbles in deionized water was 1.6 × 106, which was much lower than the prepared nanobubble water. According to previous studies [2,43], the mixing of organic solvents with pure water leads to the spontaneous formation of nanobubbles, while nanobubbles cannot form spontaneously in deionized water. the concentration of the bulk of the nanobubbles shifted to the right in the 0.05 M CaCl2 and 0.05 M Na2CO3 solutions, which indicated that their mean diameters increased. Simultaneously, in the presence of a strong electrolyte, the mono-dispersity and the total concentration of the bulk nanobubbles reduced, especially noticeable in the 0.05 M CaCl2 solution.

**Table 1.** Mean diameter, zeta potential, and concentration of nanobubbles in deionized water, in 0.05 M CaCl2 aqueous solution or 0.05 M Na2CO3 aqueous solution.


**Figure 2.** Size distribution of nanobubbles in deionized water with nanobubbles, in 0.05 M CaCl2 aqueous solution, 0.05 M Na2CO3 aqueous solution or in deionized water without nanobubbles.

#### *3.2. The Influence of Nanobubbles on the Transformation from Vaterite to Calcite*

To study the influence of nanobubbles on the transformation from vaterite to calcite, deionized water or nanobubble water was used as a solvent. The products obtained with

different reaction times in the presence/absence of nanobubbles were characterized by XRD. The XRD patterns of the obtained samples with different reaction times at 25 ◦C in nanobubble water are given in Figure 3a. In comparison with their standard JCPDS files (vaterite, 720506) and (calcite, 830577), at reaction times 0.5, 1, 2, and 3 h, a mixture of calcite and vaterite was identified. The peaks corresponding to vaterite vanished as the reaction time progressed, while the peak at 2θ = 29.46 corresponding to calcite became sharper and higher. After 4 h, pure calcite was identified. However, the XRD patterns of the obtained samples with different reaction times in deionized water in Figure 3b indicated that the peaks corresponding to calcite rapidly increased in intensity as vaterite transformed into calcite, but after 4 h, a mixture of calcite and vaterite was still identified.

**Figure 3.** XRD patterns of CaCO3 particles obtained with different reaction times in deionized water containing nanobubbles (**a**) and deionized water (**b**).

Comparing Figure 3a,b, when the deionized water containing nanobubbles was used as a solvent, the transformation from vaterite to calcite was completed in less time, indicating that nanobubbles could enhance this transformation.

According to the XRD patterns, the phase compositions calculated from Equations (1) and (2) are shown in Figure 4, which demonstrates the change in the mass fraction of calcite. Regardless of whether pure water or nanobubble water was used as the solvent, the vaterite mass fraction decreased and vanished with reaction time. After 1 h, the rate of transformation from vaterite to calcite decreased, and this might be because the nanobubbles gradually decreased with the progress of the reaction, and their effect on the transformation became weaker. After 4 h, when nanobubble water was used as the solvent, all the vaterite was converted to calcite. However, when deionized water was used as the solvent, a mixture of calcite and vaterite remained, with the mass of calcite accounting for only about 60% after 4 h, and took more than 6 h to finish this transformation from vaterite to calcite. The transformation from vaterite to calcite is a spontaneous process. The driving force for the acceleratory transformation by the present of nanobubbles was the vaterite–nanobubble coagulation. This is described in the following discussion Section 4.

To confirm whether vaterite remains in the products as the reaction progresses, FT-IR spectra of CaCO3 were obtained after 4 h in the nanobubble water and deionized water and are plotted in Figure 5. The absorption peaks at 709 cm−<sup>1</sup> indicate the presence of calcite and the peak at 748 cm−<sup>1</sup> verifies the presence of vaterite [46–48]. Figure 5 indicates the presence of calcite in deionized water containing nanobubbles (a) while the presence of a mixture of calcite and vaterite in deionized water (b), and these results were in good accordance with the XRD results shown in Figure 3. The absorption peak at 709 cm−<sup>1</sup> (inplane vibration) of calcite (b) obtained in deionized water is slightly red shifted, compared with that of sample (a). The phase transformation from vaterite to calcite might have caused the peak shift due to the slight change in the atomic distance [49].

**Figure 4.** Mass fraction of the calcite (Wcalcite/Wcalcite+vaterite) as a function of the reaction time in deionized water containing nanobubbles and deionized water.

**Figure 5.** FT-IR spectra of CaCO3 obtained after 4 h in deionized water containing nanobubbles (a) and deionized water (b).

SEM images of the fabricated samples are shown in Figure 6, with different reaction times in nanobubble water (a) to (d) and deionized water in the absence of nanobubbles (e) to (h). In deionized water containing nanobubbles, a mixture of spherical particles, i.e., crystals of vaterite [22,25,50], and cubic crystals, i.e., crystals of calcite [22,25,50], were found as shown in Figure 6a–c. It was obvious that the number of cubic crystals in deionized water containing nanobubbles increased over time with only cubic crystals present after 4 h, as shown in Figure 6d. On the other hand, in deionized water without nanobubbles, all of the products obtained at different reaction times were a mixture of spherical particles and cubic particles as shown in Figure 6e–h. The results obtained by SEM images agreed well with the results obtained by XRD (Figures 3 and 4) and FT-IR (Figure 5).

**Figure 6.** SEM images of CaCO3 obtained in deionized water containing nanobubbles at different reaction times1h(**a**), 2 h (**b**), 3 h (**c**), and 4 h (**d**) and in deionized water at different reaction times 1 h (**e**), 2 h (**f**), 3 h (**g**) and 4 h (**h**).

#### *3.3. The Influence of Concentration of Nanobubbles on the Transformation from Vaterite to Calcite*

To study the effect of nanobubbles on the polymorphs of CaCO3, deionized water with different concentrations of nanobubbles was used as a solvent. In this experiment, high-speed shearing was used to produce nanobubbles. The same motor frequency and shearing time were chosen to avoid multiple variables occurring at the same time.

The volume ratio of deionized water to nanobubbles in groups 1 to 4 is listed in Table 2. The low concentration of nanobubbles in deionized water (Groups 2 and 3) was prepared by adding deionized water containing nanobubbles obtained by high-speed shearing in different volumes to the deionized water. In groups 2 and 3, the volume ratio of pure water to nanobubbles (Vdeionized water containing nanobubble: Vdeionized water) was 1:3 and 1:1, respectively. In group 1, only deionized water was used as a solvent, while in group 4, nanobubble water without dilution was used. The total nanobubble number was from 0 to 2.96 × <sup>10</sup><sup>8</sup> bubbles/mL.


**Table 2.** The concentration of nanobubbles in deionized water in groups 1 to 4.

The XRD patterns of the obtained samples with different concentrations of nanobubbles in deionized water are shown in Figure 7. When the concentration of nanobubbles (groups 1 to 3) were low, mixtures composed of vaterite and calcite were obtained. When the concentration of nanobubbles was 2.96 × <sup>10</sup><sup>8</sup> bubbles/mL (group 4), only peaks corresponding to calcite were found.

According to the calculations using Equations (1) and (2), the mass fraction of the calcite (Wcalcite/Wcalcite+vaterite) as a function of the concentration of nanobubbles in deionized water after 4 h is shown in Figure 8. As the concentration of nanobubbles increased, the mass fraction of calcite increased linearly and reached 100% when the concentration of nanobubbles was at 2.96 × <sup>10</sup><sup>8</sup> bubbles/mL.

**Figure 7.** XRD patterns of CaCO3 particles obtained after 4 h as a function of the concentration of nanobubbles in deionized water in groups 1 to 4.

**Figure 8.** Mass fraction of the calcite (Wcalcite/Wcalcite+vaterite) as a function of the concentration of nanobubbles in deionized water after 4 h.

SEM images of CaCO3 obtained after 4 h in groups 1–4 are shown in Figure 9a–d, respectively. In Figure 9a–c, the products obtained were a mixture of spherical vaterite particles and cubic calcite particles, while in deionized water containing a high concentration of nanobubbles (2.96 × <sup>10</sup><sup>8</sup> bubbles/mL), only cubic calcite particles were observed, as shown in Figure 8. The results obtained by SEM images were quite consistent with those obtained by XRD, as shown in Figure 7. The size of calcite and vaterite particles was almost comparable in deionized water containing various concentrations of nanobubbles. The size distribution of calcite after 4 h prepared in 2.96 × 108 bubbles/mL of nanobubble water is shown in Figure S2a. The mean diameter of calcite was around 5 μm.

**Figure 9.** Mass fraction of the calcite (Wcalcite/Wcalcite+vaterite) as a function of the concentration of nanobubbles in deionized water after 4 h, group 1 (**a**), group 2 (**b**), group 3 (**c**), group 4 (**d**).

#### **4. Discussion**

*4.1. CaCO3 Particle Synthesis with Different Nanobubble Concentrations* The reaction of CaCO3 synthesis is as follows:

$$\text{Na}\_2\text{CO}\_3 + \text{CaCl}\_2 \to \text{CaCO}\_3 + 2\text{NaCl} \tag{3}$$

During this reaction, 0.05 M of Na2CO3 aqueous solution was reacted with an equal amount of CaCl2 aqueous solution of the same concentration, the synthesized CaCO3 was 0.025 M (C1). One particle diameter (d) of synthesized CaCO3 particle was around 5 μm as shown in Figure S2a. The average density of CaCO3 (*ρ*) was about 2.7 g/mL (calcite 2.71, vaterite 2.65 [51]). The total number of CaCO3 particles synthesized (NCaCO3) is as follows,

$$\text{N}\_{\text{CaCO}\_3} = \text{V/V}\_{\text{CaCO}\_3} = 6 \text{M}\_{\text{CaCO}\_3} \times \text{C}\_1 \times \text{V}\_{\text{solution}} / (\pi \rho \text{d}^3) \tag{4}$$

The concentration of CaCO3 particles (C, particles/mL) is given by Equation (5),

$$\text{C} = \text{N}\_{\text{CaCO}\_3}/\text{V}\_{\text{solution}} = 6 \text{M}\_{\text{CaCO}\_3} \times \text{C}\_1/(\pi \text{pd}^3) = 1.4 \times 10^7 \text{particles/mL} \tag{5}$$

where V is the total volume of calcium carbonate obtained, VCaCO3 is the volume of a single calcium carbonate particle, MCaCO3 is the molecular weight of calcium carbonate, Vsolution is the volume of the solution, and d is the mean diameter of the calcium carbonate particles.

As the nanobubble number was 2 to 3 × <sup>10</sup><sup>8</sup> bubbles/mL as shown in Table 2, the number of nanobubbles applied for this experiment was 10 times as many as that of CaCO3 particles.

The interaction between nanobubbles and CaCO3 particles is important for the preparation of calcite and vaterite. After the reaction shown in Equation (3), the pH equilibrated around 10. Figure 10 shows the zeta potential of air nanobubbles, calcite, and vaterite in the 0.05 M of NaCl aqueous solution containing different concentrations of Ca2+ at pH = 10. At pH = 10, the zeta potential of vaterite was always positive, while calcite and air nanobubbles were negative at concentrations less than 10−<sup>3</sup> M of Ca2+.

After the reaction described in Equation (3), the Ca2+ concentration in the solution was constant at 2.5 × <sup>10</sup>−<sup>5</sup> M after 1.5 h, as shown in Figure S4. Therefore, after 1.5 h the zeta potential of calcite and nanobubbles was negative and the zeta potential of vaterite was positive, as shown in Figure 10.

#### *4.2. The Interaction between Nanobubbles, CaCO3 Particles and Nanobubble-CaCO3 Particles*

In this section, we evaluate (a) bubble–bubble interaction, (b) vaterite–calcite interaction, (c) bubble–calcite interaction, and (d) bubble–vaterite interaction, by using the classical and extended DLVO theory in order to investigate and propose the mechanism of the transformation from vaterite to calcite in the presence of air nanobubbles.

**Figure 10.** Zeta potential of air nanobubble, calcite, and vaterite in a 0.05 M of NaCl aqueous solution containing various concentrations of Ca2+ at pH = 10.

In addition to the classical DLVO theory, i.e., the summation of the attractive van der Waals interaction energy (EA) and the repulsive electrostatic interaction energy (ER), to evaluate the interactions between different CaCO3 particles (vaterite-calcite), we used the hydrophobic interaction energy (Eh) to evaluate the stability of the nanobubbles. The total potential energy (ET) between two nanobubbles is described as follows [8,52,53]:

$$\mathbf{E\_{T}} = \mathbf{E\_{A}} + \mathbf{E\_{R}} + \mathbf{E\_{h}} \tag{6}$$

$$\frac{\mathbf{E\_T}}{\mathbf{kT}} = \frac{\mathbf{E\_A} + \mathbf{E\_R} + \mathbf{E\_h}}{\mathbf{kT}} \tag{7}$$

The total potential energy (ET) is a key indicator of the stability of nanobubbles and can evaluate whether the nanobubbles coagulate or disperse. When ET is significantly positive (e.g., >15 kT), there is a potential barrier between the two adjacent bubbles which will not coagulate, while if ET is slightly positive or negative (<15 kT), the nanobubbles can coagulate together and form a larger bubble.

EA + Eh and ER can be described by the following formulas:

$$\mathbf{E}\_{\rm A} + \mathbf{E}\_{\rm h} = -\frac{\mathbf{A} + \mathbf{K}}{6} \{ \frac{2\mathbf{a}\_1 \mathbf{a}\_2}{\mathbf{R}^2 - (\mathbf{a}\_1 + \mathbf{a}\_2)^2} + \frac{2\mathbf{a}\_1 \mathbf{a}\_2}{\mathbf{R}^2 - (\mathbf{a}\_1 - \mathbf{a}\_2)^2} + \ln\left[\mathbf{R}^2 - (\mathbf{a}\_1 + \mathbf{a}\_2)^2\right] - \ln[\mathbf{R}^2 - (\mathbf{a}\_1 - \mathbf{a}\_2)^2] \}\tag{8}$$

$$\mathbf{R} = \mathbf{a}\_1 + \mathbf{a}\_2 + \mathbf{h} \tag{9}$$

$$\mathbf{E}\_{\rm R} = \frac{\pi \varepsilon\_{\rm r} \varepsilon\_{0} \mathbf{a}\_{1} \mathbf{a}\_{2} (\boldsymbol{\psi}\_{1}^{2} + \boldsymbol{\psi}\_{2}^{2})}{\mathbf{a}\_{1} + \mathbf{a}\_{2}} \{ \frac{2 \boldsymbol{\psi}\_{1} \boldsymbol{\psi}\_{2}}{\boldsymbol{\psi}\_{1}^{2} + \boldsymbol{\psi}\_{2}^{2}} \ln \frac{1 + \exp(-\kappa \mathbf{h})}{1 - \exp(-\kappa \mathbf{h})} + \ln[1 + \exp\left(-2\kappa \mathbf{h}\right)] \} \tag{10}$$

where A is the Hamaker constant, K is the hydrophobic force constant, a1 and a2 are the radii of the nanobubbles, h is the surface distance of the nanobubbles, κ represents the Debye–Hückel parameter, and ψ<sup>1</sup> and ψ<sup>2</sup> represent the surface potential of the nanobubbles of radii a1 and a2, respectively. The thickness of the electric double layer can be represented by the Debye length (λ<sup>D</sup> = 1/κ), which can be described by the following formula:

$$\kappa = \sqrt{\frac{2\mathbf{n}\mathbf{z}^2\mathbf{e}^2}{\varepsilon\_\mathbf{r}\varepsilon\_0\mathbf{kT}}}\tag{11}$$

$$\mathbf{n} = 1000 \text{N}\_{\text{A}} \text{C} \tag{12}$$

where n is the concentration of anions or cations in the solution, z represents the valence of an ion, e represents the electron charge, ε<sup>r</sup> represents the relative dielectric constant, ε<sup>0</sup> represents the permittivity of a vacuum, T represents the temperature (K), NA represents the Avogadro number and C represents concentration of anions or cations (mol/L).

Bubble–bubble interaction: according to our previous study [8], these nanobubbles could stay stable for several weeks or even months. Even in the presence of salt, nanobubbles can also be stable for days. They were stable enough to perform zeta potential measurements. The zeta potential of nanobubbles is shown in Figure 10. When the nanobubbles existed in the reacted solution with a Ca2+ concentration less than 10−<sup>3</sup> M, the zeta potential of the nanobubble was negative, as shown in Figure 11. The mean diameter of the nanobubbles in the mixed vaterite and calcite suspension after reaction depending on the time measured by DLS after filtration through 10–15 μm filter paper is shown in Figure 11. The mean diameter of the nanobubbles increased from an initial 100–130 nm (Table 1) to 600–700 nm during 3 h of the reaction time (Figure 11). Figure 12a shows the total potential energy (ET) between the bubbles calculated by using the extended DLVO theory, and there was no potential barrier to prevent the coalescence of the bubbles, which was consistent with the bubble size enlargement as shown in Figure 11.

**Figure 11.** The mean diameter of nanobubbles in mixed vaterite and calcite suspension after reaction depending on the time measured by DLS after filtration with 10–15 μm filter paper.

Vaterite–calcite interaction: Figure 12b shows the total potential energy (ET) between the calcite and vaterite particles (6 μm mean diameter as shown in Figure S3c) calculated by using the classical DLVO theory as this was an interaction between solid particles. There was no potential barrier when considering the low zeta potential of calcite (−10 mV) and vaterite (20 mV) indicating their coagulation, (Figure 10).

Air bubble–calcite interaction: Figure 13 shows the total potential energy (ET) between the different sizes of air bubbles and calcite particles as a function of their distance by the extended DLVO theory. Figure 13a,b depict the interactions between the initial 120 nm and 700 nm nanobubbles after 1 h (shown in Figure 11) and the prepared 5 μm calcite after 1 h (shown in Figure S2b), respectively. In the both cases, there was no potential barrier between nanobubbles and calcite, indicating their hetero-coagulation.

**Figure 12.** (**a**) Total potential energy ET between bubbles (700 nm, −15 mV) as a function of the distance between two bubbles by the extended DLVO theory, at a Na+ concertation of 0.05 M and a Ca2+ concentration of 2 <sup>×</sup> <sup>10</sup>−<sup>5</sup> M; (**b**) Total potential energy ET between calcite (3 <sup>μ</sup>m, <sup>−</sup>10 mV) and vaterite (6 μm, 20 mV) as a function of the distance between two particles by the classical DLVO theory, at a Ca2+ concentration of 2.5 <sup>×</sup> <sup>10</sup>−<sup>5</sup> M.

**Figure 13.** Total potential energy ET between the different sizes of air bubbles and calcium carbonate particles as a function of the distance between a nanobubble and a calcium carbonate particle by the extended DLVO theory (K = 10−<sup>19</sup> J). Na+ concentration was 0.05 M and Ca2+ concentration was <sup>2</sup> <sup>×</sup> <sup>10</sup>−<sup>5</sup> M. (**a**) Air nanobubbles (120 nm, <sup>−</sup>15 mV) and calcite particles (5 <sup>μ</sup>m, <sup>−</sup>10 mV); (**b**) Air nanobubbles (700 nm, −15 mV) and calcite particle (5 μm, −10 mV); (**c**) Air nanobubbles (120 nm, −15 mV) and vaterite particles (6 μm, 20 mV); (**d**) Air nanobubbles (700 nm, −15 mV) and vaterite particles (6 μm, 20 mV).

Air bubble–vaterite interaction: considering the zeta potential of the vaterite particles was positive (20 mV) and air bubbles was negative, as shown in Figure 13c,d, there was no potential barrier between the nanobubbles and vaterite indicating hetero-coagulation between nanobubbles and vaterite.

To sum up, in deionized water containing nanobubbles, according to the calculations based on the extended DLVO theory, there was no potential barrier between nanobubbles, calcite or vaterite, indicating that negatively charged nanobubbles could attach and accumulate not only on the surfaces of calcite particles but also on the surfaces of vaterite particles. This supported the premise that nanobubbles could bind calcite with vaterite.

#### *4.3. The Capillary Force Model between Calcite and Vaterite*

Figure 14 shows SEM images of spherical vaterite particles attached to the cubic surface of calcite particles in deionized water (a) and nanobubble-containing solution (b) after 0.5 h. Obviously, in the nanobubble-containing solution (b), there were more spherical vaterite particles attached to the surface of cubic calcite particles than that in the deionized water without nanobubbles (a). The zeta potential of the vaterite particles was positive, while the surface of the calcite particles was negative. The vaterite particles tended to be clustered with calcite particles. Furthermore, compared to vaterite particles, calcite particles were cubic crystals, which had a higher probability of coming into contact with vaterite particles. This indicates that nanobubbles can bind calcite with vaterite, enabling vaterite and calcite to clustered together in deionized water containing nanobubbles.

**Figure 14.** SEM images of spherical vaterite particles attached onto cubic calcite particle surfaces in deionized water (**a**) and in deionized water containing nanobubbles (**b**) (after 0.5 h passed).

The capillary force model [54–57], based on the forces caused by the pressure drop across the interface at the neck and the interfacial tension force, was used to calculate the capillary force between the vaterite and calcite particles:

$$\mathbf{F}\_{\mathbf{c}} = \pi \sigma \sin \alpha [2 \sin(\alpha + \theta) + \mathrm{R} \sin \alpha (1/\mathbf{r} - 1/\mathbf{l})] \tag{13}$$

$$\mathbf{r} = \left[ \mathbf{R} (1 - \cos \alpha) + \mathbf{h} \right] / \left[ \cos(\alpha + \theta) + \cos \alpha \right] \tag{14}$$

$$\mathbf{l} = \mathbf{R}\sin\alpha - \mathbf{r}[1 - \sin(\alpha + \theta)] \tag{15}$$

where Fc is the capillary force, σ is the surface tension coefficient of the liquid, α is the fill angle, R is the radius of the colloidal probe, θ is the liquid contact angle, r and l are the principal radii of the capillary bridge, and h is the separation distance between the vaterite and calcite surfaces, as shown in Figure 15a.

In absence of nanobubbles, spherical vaterite attached to the surface of cubic calcite might become separated by the repeated collisions between solid particles and high-speed flowing solution stirred under the magnetic agitator [58]. It can be difficult for ultrafine particles to collide and coagulated. On the other hand, in presence of nanobubbles, the negatively charged nanobubbles can attach and accumulate on the surface of calcite particles and this aggregate can further coagulate with positively charged vaterite particles [58].

**Figure 15.** The geometry of vaterite and calcite particles used to calculate the capillary force between them (**a**); simplified mechanism of a nanobubble bridging the surfaces of vaterite and calcite particles (**b**); the nanobubble capillary force bridging vaterite and calcite particles as a function of their surface distance (**c**).

According to previous studies, the probability of collision between a particle and a bubble is closely related to bubble size [59]. Because of the smaller contact area and aggregation effects on ultrafine particles, nanobubbles are more easily attached to the surface of vaterite and calcite particles [13] in comparison with larger bubbles. According to the nanobubble bridging capillary force (NBCF) resulting from the long range hydrophobic attractive force [55,60–62], as they approached, a gaseous capillary bridge could be formed through the coalescence of nanobubbles. The underlying mechanism is described in Figure 15b.

The induced attractive capillary force was helpful to enhance the aggregation of vaterite and calcite particles through bubble bridging. The nanobubble bridging capillary force as a function of the surface separation distance between vaterite and calcite particles is shown in Figure 15c. Even if the separation distance was 1000 nm, Fc was still 0.001 mN/m, which would avoid separation of concentrated solid particles and enhance the effective slip lengths between solid and liquid, enhancing the aggregation between solid particles [63]. Once vaterite and calcite particles aggregated together, the capillary forces produced in the presence of the nanobubbles prevented them from separating. The calculation results agreed well with the results obtained by SEM images (Figure 14). There were more spherical vaterite particles attached to the surfaces of cubic calcite particles in deionized water in the presence of nanobubbles than that in deionized water in the absence of nanobubbles.

#### *4.4. Formation Mechanism of the Crystallization Transformation of CaCO3*

In this series of experiments reported in this article, the double decomposition method was used to produce calcium carbonate. After mixing CaCl2 and Na2CO3 solution, amorphous calcium carbonate (ACC) formed in several seconds, then ACC transformed into vaterite, and finally, all of them changed into calcite. This transformation occurred in two stages: the first step was the transformation from ACC to vaterite, and the second step was the rate-determining process from vaterite to calcite. Gas–water interfaces might act as nucleation sites, which would enhance the transformation from ACC to vaterite. According to the literature [22,23,64], the second stage of the reaction is much slower than the first stage. However, the second step directly determines the formation of calcite. Through the dissolution–reprecipitation mechanism, vaterite can be transformed into calcite [22,23,64]. For detail, the process was divided into two steps, the first step was the dissolution of vaterite, and the second step was the growth of the calcite, which was the rate-determining process. The second step was the recrystallization of calcium carbonate which was divided into the dissolution of vaterite and the growth of calcite [22,50,64]. As shown in Figure 16, Vaterite can dissolve in water and release calcium and carbonate ions, which are then reprecipitated on the surfaces of the growing calcite crystals. As shown in Figures 4 and 8, after 4 h, all the vaterite was converted to calcite in deionized water containing nanobubbles, while, in absence of nanobubbles, a mixture of calcite (60%) and vaterite (40%) remained. With the concentration of nanobubbles increasing, the mass fraction of calcite increased linearly and reached 100% when the concentration of nanobubbles was at 2.96 × 108 bubbles/mL. The results show that nanobubbles can enhance the transformation from vaterite to calcite. According to the calculation shown in Figure 13, there was no potential barrier between nanobubble, calcite and vaterite, indicating heterocoagulation between them. The number of nanobubbles applied for this experiment was 10 times as many as that of CaCO3 particles. Under the intense stirring of the magnetic agitator, nanobubbles improved the probability of collision. Due to the nanobubble bridging capillary force [63] increases with shorter separation distances between the vaterite and calcite particles, as shown in Figure 15c; thus, their aggregation can become more stable. As clearly shown in Figure 14, it was easier and more stable for vaterite and calcite to become aggregated together in the presence of nanobubbles (a) compared with the absence of nanobubbles (b). This aggregation allowed more Ca2+ to attach to the surface of calcite, allowing calcite to grow more efficiently and transform faster.

**Figure 16.** The interaction between nanobubble, vaterite, and calcite and transformation from vaterite to calcite. (Ca2+ concentration 10−<sup>4</sup> to 10−<sup>6</sup> M, pH = 10, vaterite mean diameter 5 μm, air nanobubble 120 nm.).

#### **5. Conclusions**

This was the first study using nanobubbles to aid the transformation of vaterite to calcite, and we confirmed that the transformation was accelerated in the presence of air nanobubbles. According to the calculations based on the extended DLVO theory, in the reacted solution, there was no potential barrier between the nanobubbles and both calcite or vaterite particles, indicating their coagulation. In deionized water without nanobubbles, spherical vaterite attached on the surface of cubic calcite particles of about 5 μm in size but would be dispersed by the repeated collisions between solid particles and high-speed flowing solution. However, according to the calculations of the nanobubble bridging capillary force and the SEM images of synthesized calcite particles, nanobubbles played an important role as the binder in strengthening the coagulation between calcite and vaterite. Therefore, the transformation rate increased in deionized water containing nanobubbles, compared with the aqueous solution without nanobubbles. Based on the results, this research provides a novel method for controlling the crystal morphology and reaction kinetics of CaCO3. The CO2 bubbling method (indirect carbonation) is an important way to synthesize vaterite. The size of the bubbles has an obvious effect on the vaterite content and can be controlled by adjusting the ventilation equipment. In future studies, the influence of the bubble size on crystal formation of calcium carbonate should be studied to produce pure vaterite.

**Supplementary Materials:** The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/ma15217437/s1, Figure S1: SEM images and XRD pattern of vaterite obtained; Figure S2: The schematic diagram of the machine for producing bulk nanobubbles (a), and the size of the stainless-steel blade and beaker (b); Figure S3: Size distribution of calcite after 4 h (a), the size distribution of calcite after 1 h, and size distribution of vaterite after 1 h prepared in 2.96 × 10<sup>8</sup> particles/mL of nanobubble water; Figure S4: The concentration of Ca2+ depending on the time in the solution filtrated at 0.45 μm after the reaction to prepare CaCO3;

**Author Contributions:** Y.W.: Data curation, Investigation, Writing original draft. M.H.: Investigation. C.H.: Review and editing. K.W.: Review and editing. N.T.H.N.: Writing—Review and editing. S.L.: Review and editing. G.D.: Review and editing. A.O.: Review and editing. T.F.: Writing—Review and editing. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Natural Science Foundation of China, grant number NSFC 21976039.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The authors confirm that the data supporting the findings of this study are available within the article.

**Acknowledgments:** We appreciate the special funding support by "Guangxi Bagui scholars".

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

#### **References**


## *Article* **Simulation of Electrical and Thermal Properties of Granite under the Application of Electrical Pulses Using Equivalent Circuit Models**

**Kyosuke Fukushima 1, Mahmudul Kabir 1,\*, Kensuke Kanda 1, Naoko Obara 1, Mayuko Fukuyama <sup>2</sup> and Akira Otsuki 3,4**


**Abstract:** Since energy efficiency in comminution of ores is as small as 1% using a mechanical crushing process, it is highly demanded to improve its efficiency. Using electrical impulses to selectively liberate valuable minerals from ores can be a solution of this problem. In this work, we developed a simulation method using equivalent circuits of granite to better understand the crushing process with high-voltage (HV) electrical pulses. From our simulation works, we calculated the electric field distributions in granite when an electrical pulse was applied. We also calculated other associated electrical phenomena such as produced heat and temperature changes from the simulation results. A decrease in the electric field was observed in the plagioclase with high electrical conductivity and void space. This suggests that the void volume in each mineral is important in calculating the electrical properties. Our equivalent circuit models considering both the electrical conductivity and dielectric constant of a granite can more accurately represent the electrical properties of granite under HV electric pulse application. These results will help us better understand the liberation of minerals from granite by electric pulse application.

**Keywords:** electric field; temperature distribution; conductivity; dielectric constant; liberation; hardrock; mineral distribution

#### **1. Introduction**

Minerals are familiar to us in our daily lives and are used in various fields (e.g., ceramics, pharmaceuticals, cosmetics) [1]. In order to extract minerals from rocks, a combination of mechanical comminution methods of crushing and grinding (e.g., ball mill) based on compression, impact and shearing mechanisms are generally used [2]. The energy efficiency of these methods is as low as 1%. Especially for hard rocks, such as granite, the mechanical grinding is inefficient due to tool wear and the increase in cost [3]. In addition, the energy used for comminution is very large, accounting for 2–4% of the world energy consumption [4,5]. Therefore, a new method for size reduction and selectively liberating minerals from rocks using high-voltage electric pulses has been developed [2]. Electric pulse liberation of minerals has been reported as a solution for low energy efficiency in comminution [6–8].

**Citation:** Fukushima, K.; Kabir, M.; Kanda, K.; Obara, N.; Fukuyama, M.; Otsuki, A. Simulation of Electrical and Thermal Properties of Granite under the Application of Electrical Pulses Using Equivalent Circuit Models. *Materials* **2022**, *15*, 1039. https://doi.org/10.3390/ ma15031039

Academic Editor: Georgios C. Psarras

Received: 31 December 2021 Accepted: 26 January 2022 Published: 28 January 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/).

There are two types of electrical disintegration methods, and they are (1) Electro Hydraulic Disintegration (EHD) and (2) Electrical Disintegration (ED). Figure 1 shows a schematic of the EHD and ED methods. For (1) EHD, an upper/high voltage (HV) electrode and a ground electrode are placed above the object in a liquid and a shock wave is generated for the selective liberation of minerals (Figure 1a). (2) The ED method is also applied in a liquid but the upper/HV electrode is placed directly on or very close to the object to be disintegrated (Figure 1b). The ground electrode is placed below the object and when the electric pulse is applied, a large electric current is generated and flows through the object and thus leads to its disintegration. Water, due to its low electrical conductivity, is usually used as a liquid. For the EHD method, the voltage is applied only to the rock because the breakdown voltage of the rock is lower than the electric field strength of the water when the rise time of the electric pulse is less than 500 ns [9]. In the EHD method, lightning impulses are applied to two electrodes in close proximity, and the compressive force generated by the shock wave is used for crushing. In the ED method, the upper/HV electrode is placed close to the sample, and an instantaneous large current is applied to the inside of the sample to crush it [9]. The ED method is used in this study because of its potential for selective grinding and its ability to weaken the rock and reduce the energy consumption for grinding hard rocks [10,11].

**Figure 1.** Model diagrams of different electric pulse crushing [4]. (**a**) EHD method; (**b**) ED method.

The ED method has been applied to the separation of precious minerals such as diamonds and emeralds [11]. In addition, the effects of electric pulse application on rocks have been studied based either on the electric conductivity or dielectric constant of rocks experimentally or simulation [11–13]. Walsh et al. (2020) performed simulation of hard rocks under HV pulse application considering the microstructures (i.e., distribution of minerals) of the rocks and they described about the change in minerals of the rocks under HV application. Li et al. (2018) showed how the *I*-*V* properties of granite changed under HV applications with different types of electrode spacing. However, the effects of electric pulses on rocks are still largely empirical, and our work is the very first study that considered both conductivity and dielectric constant as far as the authors are aware of.

The purpose of this study is to obtain better knowledge of granite comminution and selective liberation under high voltage pulses. Granite was selected as an example of hard rock. For this purpose, we developed an equivalent circuit model of granite based on both conductivity and dielectric constant, considering the minerals and voids in the granite, and applied a high electric field to the granite model to simulate electrical characteristics under HV applications. Specifically, our equivalent circuit model was created using a C# program, and the electric field distribution in granite was calculated using the circuit simulation software LTspice. Since rocks consist mostly of dielectric materials, our simulation works can be significantly helpful to understand the behavior of rocks (i.e., granite) under HV electric pulse, as the circuit models were proposed considering both electric and dielectric properties of granite.

#### **2. Simulation Methods**

#### *2.1. Simulation Model*

The simulation was performed for a well-known hard rock, i.e., granite composed of quartz, plagioclase, k-feldspar, and biotite [12]. In this study, a granite specimen with a cube of 1 mm on each side was considered, assuming that the entire rock is divided into an arbitrary number of smaller cubes to create an equivalent circuit model. In general, the dimension of a rock (i.e., granite in this work) is larger than that of our simulation works. We assumed the granite as a small cube of 1-mm size on each side in order to avoid the complexity of the work. As the electrical properties such as conductivity and relative permittivity do not vary with their sizes, it is quite reasonable to consider a smaller size for simulation works. By comparing larger and smaller sizes of the rocks, the results will show similar trends and this point was reconfirmed from our simulation results, similar to what was reported in the literature [12]. To obtain accurate results, we assumed 1000 parts of smaller cubes which make the cube of 1-mm edge length of granite, as shown in Figure 2a. The electrical properties of minerals can be represented by a parallel circuit of resistors and capacitors (Figure 2b). The resistance of minerals in the equivalent circuit model is calculated from the following equation:

$$R = \frac{1}{\sigma} \frac{L}{S} \tag{1}$$

where *R* is the resistance, *σ* is the electrical conductivity, *L* and *S* represent the length of the segmented mineral (0.1 mm) and the surface area of the segmented mineral (0.01 mm2), respectively. Capacitance is expressed by the following equation:

$$\mathbf{C} = \varepsilon \frac{\mathbf{S}}{L} \tag{2}$$

here, *C* is capacitance and *ε* is dielectric constant. The equivalent circuit was created by assigning the minerals contained in the granite to the divided cubes, calculating the resistance and capacitor values for each mineral, and constructing the circuit. As an example, Figure 2 shows the equivalent circuit model for a 10 × 10 × 10 partitioned cube. The simulation work was performed for a needle-flat electrode as shown in Figure 2a. From the knowledge of electric circuits, it can be understood that the electric field strength is mainly applied in the vertical direction. The electric field strength in the horizontal direction was neglected because the values were less than 3% compared to the ones in vertical direction. Thus, in this paper, we have shown the electric field strengths in the vertical direction only. The potential difference for each circuit element was calculated by dividing the values by the distance (i.e., 0.1 mm). Then, the electric field was calculated for each smaller cube and the data were used to draw the 3D mapping of the electric field distribution with MATLAB (R2016a, The MathWorks, Inc., Natick, MA, USA). Other simulation results were also used to draw 3D images by MATLAB.

#### *2.2. Composition of Granite*

The mineral composition of a typical granite is referred to the information available in [9], as shown in Table 1. The distribution of each mineral in the granite was determined by a random placement program based on the compositional ratio using the programming software C#.


**Table 1.** Mineral composition of a typical granite [12].

**Figure 2.** Concept of equivalent circuit model. (**a**) Granite; (**b**) equivalent circuit model for a small divided cube of granite.

#### *2.3. Electrical and Thermal Properties of Minerals*

The current *J* per unit area flowing through a rock is expressed by Ohm's law:

$$\mathbf{J} = \sigma \mathbf{E} + \frac{\partial \mathbf{D}}{\partial t} = \sigma \mathbf{E} + \varepsilon \frac{\partial \mathbf{E}}{\partial t} \tag{3}$$

where *E* is the electric field, and *ε* is the dielectric constant [14]. From Ohm's law, the current distribution in a rock is affected by the electrical conductivity and dielectric constant. The relationship between temperature and electrical conductivity of minerals (*σs*) is empirically expressed by the following equation [12]:

$$
\sigma\_s = A \exp\left(\frac{-B}{k\_b T}\right) \tag{4}
$$

 where *A* and *B* are constants, *T* is the absolute temperature, and *kb* is Boltzmann's constant (8.618 × <sup>10</sup>−<sup>5</sup> eV/K). The constants *<sup>A</sup>* and *<sup>B</sup>* of the granite rocks are shown in Table <sup>2</sup> [12], and the electrical conductivity of the minerals as a function of temperature is shown in Figure 3. The values of conductivity in Figure 3 were calculated for each single mineral using Equation (4). It is worth mentioning that the conductivity values of minerals (i.e., Figure 3) are surely different from the conductivity values of a mix solution (i.e., for NaCl solution of Equation (7)). At low temperature, the electrical conductivity of minerals is mainly affected by the water content in the fractures and pores. Sinmyo and Keppler reported the following empirical equation for NaCl solutions [15]:

$$
\log\_{10} \sigma\_{\text{Brin}} = -1.7060 \, - \, \frac{93.78}{T} + 0.8075 \log\_{10} c + 3.0781 \log\_{10} \rho + \log\_{10} (\lambda\_0(\rho, T)) \tag{5}
$$

where *c* is the mass percent concentration of NaCl, *ρ* is the mineral density, and *λ*<sup>0</sup> is the molar conductivity of NaCl in water at unlimited dilution. Here, *λ*<sup>0</sup> is approximated by the following equation:

$$
\lambda\_0 = 1573 \ -1212 \rho + \frac{537,062}{T} - \frac{208,122,721}{T^2} \tag{6}
$$

Equation (5) can only be applied to temperatures above 100 ◦C, but the following equation can be used to approximate the electrical conductivity below 150 ◦C [16]:

$$
\sigma\_f = \sigma\_{ref} \frac{T - 251.65}{T\_{ref} - 251.65} \tag{7}
$$

where *σref* is the electrical conductivity at the reference temperature, and the unit of *σref* is K. In this simulation, the *σref* , *T* and *c* was set to 150 ◦C (423.15 K), 20 ◦C (293.15 K), and 3 wt%, respectively [12,15]. Note that the temperature variation occurs instantaneously and we used Equation (17) to calculate the temperature of the granite under HV impulse application.

**Table 2.** Values of the empirical parameters [12].

**Figure 3.** Calculated values of conductivity of minerals versus temperature.

The electrical conductivity for a composite mixture of solid and fluid can be expressed using a modified Archie's law proposed by Glover [17]:

$$
\sigma\_{\rm mix} = \sigma\_s \left(1 - \mathcal{Q}\_f\right)^p + \sigma\_f \mathcal{Q}\_f^m \tag{8}
$$

where *σ<sup>s</sup>* is the electrical conductivity of a solid, *σ<sup>f</sup>* is the electrical conductivity of a fluid, ∅*<sup>f</sup>* is the volume fraction of a fluid, *<sup>m</sup>* is Archie's exponent, and *<sup>p</sup>* is determined by Equation (9).

*p* = log <sup>1</sup> <sup>−</sup> <sup>∅</sup>*<sup>f</sup> m* log <sup>1</sup> <sup>−</sup> <sup>∅</sup>*<sup>f</sup>* (9)

The value of Archie's index (*m* = 1.5) was used in this calculation and obtained from [12].

As shown in Equation (1), the value of dielectric constant is important in considering the current distribution within the minerals composing a rock. The volume of the fluid (considered equivalent to the void volume) and the dielectric constant of each mineral are shown in Table 3 [12,13].

**Table 3.** Values of void volume \* and permittivity of minerals in granite [12,13].


\* Volume filled with water during ED method application.

#### *2.4. Applied Voltage*

In the case of ED crushing, the electrodes are placed directly on the rock and the voltage is applied, and thus, the mineral crushing is theoretically possible even without using lightning impulses. However, in practice, the application of high voltage for a long period needs a large amount of electricity and lowers the energy efficiency. Due to this point, it is reasonable to use lightning impulses for mineral crushing. A lightning impulse (0.3 μs/50 μs) was used as the input voltage applied to the granite in this work. The peak value of impulse voltage was at the front time (i.e., 0.3 μs = 300 ns), and 50% of the peak voltage was at the tail time (i.e., 50 μs). In this work, the lightning impulse was simulated and applied to the granite by the circuit simulator LTspice (Figure 4). Fujita et al. (1999) explained that the voltage was applied only to the rock when the rise time of the applied voltage was less than 500 ns in the ED method [9]. Therefore, in this study, the applied voltages were: (a) 100% rise, the maximum voltage 20 kV at a rise time of 300 ns, and (b) 50% rise, 50% of the maximum voltage 10 kV at 38 ns. The electric current *I* was also calculated and plotted within the same graph of applied voltage (i.e., Figure 4). The maximum current (i.e., 2.55 mA) was observed at 2.34 ps. At the maximum value of the lightning impulse voltage (i.e., 20 kV), the value of the current was 1.53 mA at 0.3 μs. After 0.3 μs, the value of the current gradually decreased to 0.75 mA at 50 μs.

**Figure 4.** Applied voltage used in this simulation work and electric current.

#### *2.5. Temperature in the Rock*

Assuming a rock is regarded as a dielectric material, we will discuss the electric phenomena when an AC input voltage *v* with angular frequency *ω* is applied on the rock. Figure 5 shows the equivalent circuit for a typical dielectric material with AC input voltage. The AC analysis for this circuit can be explained with the use of vector analysis. The input voltage was chosen as the reference vector. The effective value of voltage *v* is *V* and the

vector expression we used is the symbol of . *V*. For the circuit shown in Figure 5a, the total electric current . *Ip* can be expressed by a vector sum of . *IC* and . *IR*. Here, . *IC* represents the displacement current and . *IR* represents the loss current of the granite. . *IC* and . *IR* can be explained by Equations (10) and (11), respectively.

$$
\dot{I}\_{\mathbb{C}} = \dot{\jmath}\omega \mathbf{C} \dot{V} \tag{10}
$$

$$
\dot{I}\_R = \frac{1}{R}\dot{V} \tag{11}
$$

Here, *<sup>ω</sup>* is angular frequency of the AC voltage *<sup>v</sup>*. The sum of the electric current, i.e., . *IP* can be found by adding Equations (10) and (11).

$$
\dot{I}\_P = \dot{I}\_R + \dot{I}\_\mathbb{C} = \left(\frac{1}{R} + j\omega\mathbb{C}\right)\dot{V} \tag{12}
$$

The amount of heat per unit time generated by the applied voltage and current flowing through the granite (i.e., electric power) can be calculated by the following equation:

$$\frac{\partial q}{\partial t} = VI\_P \cos \theta = VI\_R = \frac{1}{R}V^2 \tag{13}$$

where *q* is the heat generated, *V* is the potential difference and *θ* is the phase angle between the voltage and electric current . *IP*. Assuming the phase angle between . *IC* and . *IP* is *δ*, the relationship in between *θ* and *δ* can be expressed by the following equation,

$$
\theta + \delta = \frac{\pi}{2} \tag{14}
$$

From Figure 5b, the dielectric loss tan*δ* can be expressed by the following equation.

$$\tan \delta = \frac{\dot{I\_R}}{\dot{I\_C}} = \frac{1}{R \omega \mathcal{C}} \tag{15}$$

The values of resistance *R* and capacitance *C* were calculated for each mineral as shown in Table 4. The tan*δ* was also calculated at 20 kHz as shown in Table 4. From Equation (15), it is understood that the dielectric loss (i.e., tan*δ*) changes with frequency. The value of tan*δ* was 175 for quartz (see Table 4), but it can be changed to 1.17 at the frequency of 3 MHz. The temperature change caused by this heat can be expressed using a thermal resistance defined by the following equation:

$$R\_j = \frac{1}{k} \frac{L}{S} \tag{16}$$

where *Rj* is the thermal resistance and *k* is the thermal conductivity. The relationship between temperature change and heat is expressed by the following equation:

$$\frac{\partial q}{\partial t} = \frac{T\_1 - T\_2}{R\_j} \tag{17}$$

where *T*<sup>1</sup> and *T*<sup>2</sup> are the temperatures before and after the application of high voltage (i.e., heat generation) to the granite, respectively. The thermal conductivity of each mineral is summarized in Table 5 [18]. Equation (17) was used to calculate the temperature of granite under the HV application to be discussed later.

**Figure 5.** Equivalent circuit model for typical dielectric material. (**a**) Equivalent circuit; (**b**) vector analysis of the circuit.

**Table 4.** Calculated values of resistance *R*, capacitance *C* and tan*δ* of minerals in granite.


δ

⋅

**Table 5.** Heat characteristics of minerals in granite [18].


#### **3. Results and Discussion**

#### *3.1. Distribution of Minerals in Granite*

Using a program developed in C#, the distribution of minerals was determined, based on the mineral composition of the granite shown in Table 1. The layout information of the granite rock is shown in Figure 6a,b. Since the placement information is three-dimensional, we displayed in a cross-sectional view where the electrode is assumed to be placed, as shown in Figure 6b that is also a reference to Figure 6c–h in terms of mineral and electrode placement. The mineral distribution was chosen randomly meeting the conditions of Tables 1 and 3. Some works [19,20] reported mineral distributions arranged with Voronoi networks in order to describe any random distribution of particles, minerals, etc. Though our simulation works used a 5 times larger void, comparing to their works, the proportion of void and minerals was maintained as described in Table 3.

**Figure 6.** Simulation result when time changes. The red boxes in (**b**,**d**,**f**) are added for the visual aid. (**a**) Mineral distribution; (**b**) cross-section surface viewed from y-z direction; (**c**) electric field (50% rise of lightning impulse voltage); (**d**) electric field (100% rise of lightning impulse voltage); (**e**) heat (50% rise of lightning impulse voltage); (**f**) heat (100% rise of lightning impulse voltage); (**g**) variation of temperature (50% rise of lightning impulse voltage); (**h**) variation of temperature (100% rise of lightning impulse voltage).

#### *3.2. Simulation Results*

The rock was divided into 1000 small blocks (10 × 10 × 10), and minerals were placed in the blocks based on the mineral composition of the granite rock (Table 1). An equivalent circuit model was created based on the mineral distribution information (Figure 6b), and was used to investigate the electric field distribution in the granite when an electric pulse was applied to the granite. The applied voltage was set to reach the maximum voltage 20 kV at around 300 ns, as shown in Figure 4. The heat distribution in the granite was calculated from the electric field and the value of resistors. The temperature change was calculated from the heat distribution and the values of thermal resistance in the granite. The results of each calculated parameter for 20 kV (100% rise) and 10 kV (50% rise) are summarized in Figure 6. Table 6 shows the maximum and minimum values for each result. All of the maximum values in Table 6 were found in the small block in contact with the upper/HV electrode, and the values decreased as they moved away from the upper/HV electrode, with the minimum value found in the area in contact with the ground electrode.


Figure 6c shows the electric field distribution at the moment when 50% rise voltage was applied, and Figure 6d shows the electric field distribution at the moment when 100% rise voltage was applied. As mentioned above, the electric field distribution was high near the upper/HV electrode. In addition, the electric field distribution was found lower inside the plagioclase compared to other minerals. In Figure 6d, for example, plagioclase

placed in the area defined by the y-z coordinates (4,5), (4,6) (5,6), (5,5) showed a lower electric field of plagioclase compared to their surroundings. Plagioclase has a higher void volume than other minerals (1.8 > 0.9%) and contains more water in the ED method (see Table 3) and that can be the reason behind the lower electric field distribution in plagioclase. In this study, the initial temperature of the minerals was set to 20 ◦C, so the electrical conductivity of the mineral-water composite was highly dependent on the water content. In other words, as the amount of water increases, the electrical conductivity increases (Equation (8)). Since electrical conductivity and resistivity have an inverse relationship, the resistance decreases, and thus, the electric field decreases according to Equation (3). This might be the reason for the low electric field inside plagioclase. When the resistance of the minerals adjacent to the plagioclase is lower than that of plagioclase, the voltage of the plagioclase with lower resistance is inevitably lower than that of the surrounding area because the current flowing in the circuit branches in the vertical direction is the same. Since the electric field is calculated by dividing the potential difference between two adjacent contacts by the distance between the two contacts, the electric field of the plagioclase naturally becomes lower than the surrounding area, while a higher electric field appears in the surrounding area. In comparing 50% between 100% rise of lightning impulse voltages (Figure 6c 50% rise vs. Figure 6d 100% rise), the latter case shows a clearer difference in the electric field distribution possibly due to the change of the electric field distribution in the granite depending on the magnitude of the applied voltage. It can be due to the difference in the dielectric breakdown electric field between water and rock (minerals) and its threshold duration of the HV wavefront is about 500 ns [9]. In other words, below 500 ns, the current flows inside the minerals, otherwise it flows through the water and not through the minerals.

Figure 6e,f shows the heat distribution at the moment of applying 50% rise voltage and 100% rise voltage, respectively. We can see the change in the heat generated in the granite from Figure 6e,f calculated from Equation (13). The heat generated in the upper part of the granite was higher than that of the lower part and was similar to the electric field distribution. There was almost no heat generated in the bottom parts in granite from our simulations, especially for the 50% rise of lightning impulse voltage (Figure 6e). As shown in Equation (13), the amount of heat generated is determined by the electric field distribution and the value of resistance, and since it is known from Figure 6c,d that the resistance of plagioclase is small due to its high electric conductivity, the electric field is low. Figure 6g,h shows the temperature change calculated at 50% rise and 100% rise of lightning impulse voltage. The maximum temperature change was about 3182 K near the upper/HV electrode, but the temperature change near to the ground electrode was almost zero. In particular, the plagioclase present in the area defined by the y-z coordinates (4,5), (4,7) (6,7), (6,5) showed less heat generation than their surrounding areas which showed a relatively high heat generation compared to plagioclase. The results showed that the temperature change differs with minerals and distance from the upper/HV electrode. By comparing Figure 6e (50% rise) with Figure 6f (100% rise), the temperature distribution is more heterogeneous in the latter case and affected by the different mineral placement that can deviate the electric field distribution (Figure 6d). Again, this temperature distribution difference can be the difference in the dielectric breakdown electric field as in the case of the electric field distribution.

As seen in Figure 6g,h, the temperature change is larger in the vicinity of the upper/HV electrode, indicating that the temperature change was occurring more likely by the heat generated due to the electric field and the resistance value of each mineral rather than by the heat resistance. The results obtained in this study reproduced similar trends as reported in other research works where the calculation was performed by considering only the electrical conductivity of minerals [12,13].

Walsh et al. (2013) used Voronoi tessellation to simulate granite in order to understand the thermal spallation. Their works using the Euler–Godunov model showed that the heat distribution differs in the boundary lines of the Voronoi tessellation (i.e., minerals). Our simulation works also showed that in spite of electric field distribution differences in the granite from their work, the heat variation was very small (see Figure 6f) [19].

Walsh et al. (2020) described simulation results with a 10-cm wide granite model with thickness of 2 cm. We used the same mineral composition in the granite as they did. They performed electric pulse crushing simulation using the ED method with needle electrodes placed at the center of the granite model for both the upper/HV and ground electrodes in 2D [12]. In addition, a high electric field was observed around the electrodes. In comparison with our results, they also observed a decrease in the electric field at the plagioclase and high electric field distribution near the electrodes [12]. In our work, we used both electric and dielectric parameters to understand the behavior of granite under HV electric pulse application. From our work, as we described before, the electric properties affected the electric distribution and heat generated in the minerals. At the same time, by considering dielectric constant in our work, the reason of the generated heat difference in different minerals was clear. On the other hand, Walsh et al. (2020) used 10 times larger input voltage compared to our work in order to simulate the thermal differences. If we consider the economic side of the mineral liberation by electric pulses, our method can be used to understand electrical and thermal properties of minerals under a lower input voltage application. Within this framework, further simulation works are to be carried out to confirm the heat difference under larger HV applications.

An application example of the ED method was described in the literature where it was applied on several minerals (i.e., quartz, pyrite, calcite, albite) which were embedded in cement and their selective liberation behavior was observed [19]. The size of the cemented minerals was 10 mm × 10 mm × 8 mm and the HV electrode was a needle electrode and flat electrode was placed as a grounding electrode. Channels were formed near the contact point of the electrodes, and crushing occurred when the electric pulse was applied on the cemented minerals [21]. The formation of channels and their connections which led to crushing occurred near the needle of the electrode and was due to the high electric field generated by the electric pulse. For the sample of quartz and cement, crushing occurred along the boundary between quartz and cement. The electric field distribution changed at the boundary of cement and quartz due to the difference in their electric conductivities, and this phenomenon is similar to the relationship between plagioclase and other minerals found in our study (Figure 6d,f). The electric field data obtained in this study are consistent with the general electric-pulse comminution theory, which states a difference in electric field occurs at boundaries with large differences in electric conductivity and dielectric constant, causing rock crushing due to rapid expansion of the discharge channels and associated heat generation [2]. The above comparison with the literature indicated that our simulation results showed the same/similar trend to the previous simulation and experiments of applying ED method on rocks, and thus, our methods and results obtained in this study are considered valid.

#### **4. Conclusions**

In this study, we investigated the electrical and thermal characteristics of granite under the application of electric pulses. An equivalent circuit model of a granite in 3D (1 mm × 1 mm × 1 mm) was developed. The electric field in the circuit was then calculated considering an electric pulse application. The heat and temperature changes generated in the minerals with the granite were calculated from that result. Until now, existing studies used either electrical conductivity or dielectric constant, but there has been no model that considers both. In this study, however, we used both parameters in the equivalent circuit model, calculated the values of the electric field distribution, and then, the heat distribution, and created the temperature change distribution from the heat distribution. This has not been performed before, and can lead to more detailed data under relatively low applied voltage by considering both parameters; so, the findings of this study can be useful for energy efficient rock crushing. The results of this study are summarized as follows:


**Author Contributions:** Conceptualization, M.K. and A.O.; methodology, K.F., M.K., K.K., N.O., M.F., A.O.; software, K.F., M.K., K.K.; validation, K.F., M.K., K.K. and A.O.; formal analysis, K.F., M.K., K.K.; investigation, K.F., M.K., K.K.; resources, M.K. and A.O.; data curation, K.F., M.K., K.K. and A.O.; writing—original draft preparation, K.F., M.K.; writing—review and editing, K.F., M.K., A.O.; visualization, K.F., M.K., K.K.; supervision, M.K. and A.O.; project administration, M.K. and A.O.; funding acquisition, M.K. and A.O. All authors have read and agreed to the published version of the manuscript.

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

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

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

#### **References**


## *Article* **Equivalent Circuit Models: An Effective Tool to Simulate Electric/Dielectric Properties of Ores—An Example Using Granite**

**Kyosuke Fukushima 1, Mahmudul Kabir 1,\*, Kensuke Kanda 1, Naoko Obara 1, Mayuko Fukuyama <sup>2</sup> and Akira Otsuki 3,4,5**


**Abstract:** The equivalent circuit model is widely used in high-voltage (HV) engineering to simulate the behavior of HV applications for insulation/dielectric materials. In this study, equivalent circuit models were prepared in order to represent the electric and dielectric properties of minerals and voids in a granite rock sample. The HV electric-pulse application shows a good possibility of achieving a high energy efficiency with the size reduction and selective liberation of minerals from rocks. The electric and dielectric properties were first measured, and the mineral compositions were also determined by using a micro-X-ray fluorescence spectrometer. Ten patterns of equivalent circuit models were then prepared after considering the mineral distribution in granite. Hard rocks, as well as minerals, are dielectric materials that can be represented as resistors and capacitors in parallel connections. The values of the electric circuit parameters were determined from the known electric and dielectric parameters of the minerals in granite. The average calculated data of the electric properties of granite agreed with the measured data. The conductivity values were 53.5 pS/m (measurement) and 36.2 pS/m (simulation) in this work. Although there were some differences between the measured and calculated data of dielectric loss (*tanδ*), their trend as a function of frequency agreed. Even though our study specifically dealt with granite, the developed equivalent circuit model can be applied to any other rock.

**Keywords:** conductivity; dielectric constant; hard rock; mineral distribution; voltage-dependent resistance (VDR)

#### **1. Introduction**

The modeling of rocks is one of the effective tools available to understand the physical phenomena in rocks under different applications, including the electric-pulse liberation of minerals. The electric-pulse liberation of minerals is an emerging comminution method that can be a good solution for low-energy efficiency in size reduction and selective liberation [1,2]. As rocks are dielectric materials [3–6], it is necessary to understand their electric and dielectric properties to evaluate and better understand the physics and behaviors of minerals upon the application of a HV electric pulse on them. Since the 1990s, computer simulations have become an important tool to model the research due to the unprecedented development of

**Citation:** Fukushima, K.; Kabir, M.; Kanda, K.; Obara, N.; Fukuyama, M.; Otsuki, A. Equivalent Circuit Models: An Effective Tool to Simulate Electric/Dielectric Properties of Ores—An Example Using Granite. *Materials* **2022**, *15*, 4549. https:// doi.org/10.3390/ma15134549

Academic Editor: Saeed Chehreh Chelgani

Received: 8 April 2022 Accepted: 24 June 2022 Published: 28 June 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 data processing power of personal computers. There are many works related to electricpulse liberation being carried out using computer simulations in order to understand the HV electric-pulse liberation of rocks [7].

Andres et al. (2001) compared the electrical pulse application with the conventional mechanical comminution of oxide ores containing hematite or platinum-group metals [1]. They also showed some simulation works of the electric-field distribution in ores using the electromagnetic theory of dielectric materials.

Some other works were performed to understand the electric-pulse comminution by using the COMSOL Multiphysics software package [8,9]. Seyed et al. (2015) worked on phosphate ore that was under electric-pulse application, which showed that the electrical field was dependent on the electrical properties of minerals, particle size and the location of conductive minerals [8]. Li et al. (2018) used the COMSOL Multiphysics software to understand the behavior of electro-pulse boring in granite under a HV application, considering the composition of granite, electrode spacing and electrode shape. They found that HV boring is affected by the composition of granite and its electric properties. However, they used equivalent circuit models of an electric pulse of a HV source only [9].

Zuo et al. (2015) discussed high-voltage pulse (HVP) breakage models of ores with three breakage indices (i.e., body breakage probability, body reduction evaluation index, body breakage product pre-weakening degree) [10]. They considered the breakage results of three ore samples (i.e., gold–copper ore, iron oxide copper–gold ore (IOCG), and hematite ore) and fitted the measured data statistically. Their work linked these parameters with the mass-specific impact energy, and even though the ores were different, the HVP model showed similar results with measured data [10].

Walsh and his groups conducted some simulation works related to the HV pulse liberation of ores [11,12], using Voronoi tessellation in order to simulate the mineral distribution in the granite rocks. Their works considered the minerals of granite to understand their behaviors under HV electric applications. However, these simulation works did not consider both the electric and dielectric properties of each mineral consisting of rock.

In HV engineering, equivalent circuit models that represent insulation and dielectric materials are used to simulate and understand the behavior of HV in testing objects [9,13–17]. Zuo et al. (2020) modeled insulation cable with equivalent circuit models by considering placing the cable into small fractions of insulation (dielectric) material and discussed the HV direct current influence on insulation cables [14]. Kabir et al. (2011) [18] made equivalent circuit models to understand the fine ceramics (ZnO varistors) by considering hundreds of ZnO grains of an average grain size of 1–2 μm to understand the nonlinear *I*–*V* (current–voltage) properties of ZnO varistors. Ono et al. (2009) [17] calculated the electric-field distribution in composite materials using equivalent circuit models under HV applications. Their works discussed different types of electric pulses and electrode gaps to understand HV behavior. With the above literature, the importance of considering each mineral to understand mineral liberation under a HV application is clearly identified. Additionally, it is worth noting that the equivalent circuit model analysis has a strong potential to understand the electric pulse comminution [18].

The equivalent circuit model is used in HV applications, as discussed before. The model can simulate the behavior of the subject under a HV application effectively, but only if the equivalent circuit can imitate the subject's electrical/dielectric properties accurately. Thus, it is necessary to confirm the reproduction ability of the electrical/dielectric properties of a subject. However, the dielectric properties of composite materials are difficult to simulate by equivalent circuit models, as they are frequency-dependent and the dielectric values are fixed in an equivalent circuit model. We will discuss this regarding the simulation results.

We discussed the HV application on granite in the other paper [18], where an equivalent circuit was used to simulate the behavior of granite under a HV impulse application. As the results proved the effectiveness of our equivalent circuit model in electric-pulse liberation, in the current work, we aimed to further advance this model by considering

voids' (pores in minerals as well as rocks) behavior under a HV impulse on the rocks by using the measurement data of electrical properties of the dielectric breakdown of air. Our previous paper also dealt with voids (pores in minerals), but we assumed those voids were filled with liquid, as the electrical disintegration (ED) method applies HV impulses in liquid. Thus, the values of electric resistance that include voids, did not represent the dielectricbreakdown properties of voids (pores) in minerals; rather, it showed the electric resistance of minerals including voids, as they were assumed to be filled with liquid. Generally, the electric and dielectric properties of any dielectric and insulation materials are measured in the HV branch in air [13]. In this paper, we considered all the voids (pores in minerals), and their equivalent circuits were placed along with the circuits of minerals in the equivalent circuit model of granite. Again, we considered both the electric and dielectric properties of each type of mineral of granite (i.e., quartz, plagioclase, K-feldspar, and biotite) and made an equivalent circuit to simulate the granite's electric and dielectric properties. First, the electric and dielectric properties of the granite sample and the modal mineral composition of the sample were measured. Considering the granite ore is composed of the minerals and the voids, we calculated their values of electric resistance and capacitance and then created the circuit models to simulate the granite rock. For the simulation works, we prepared 10 types of mineral distributions randomly, and an equivalent circuit model was prepared for each mineral distribution. The simulation results of our works showed good agreement with the measured data of *I*–*V* properties and dielectric properties.

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

#### *2.1. Sample*

Our method of using an equivalent circuit to simulate rocks can be applied to any type of rock. We selected granite as an example of hard rock for this study. Our sample granite was bought from Kenis, Japan, and the sample was collected from Akaiwa city of Okayama prefecture, Japan. The sample granite was thinly sliced (with a thickness of about 1.1 mm) and its surface was well polished for measurements, including the ones for the elemental compositions as well as the electrical and dielectric properties. Figure 1a shows a snapshot of the granite sample. The elemental mapping was obtained by using a micro-X-ray fluorescence spectrometer (μXRF, M4 Tornedo Plus; Bruker, Billerica, MA, USA) and the mineral composition was then calculated by stoichiometry. Based on the imaging and X-ray analysis, whose results are shown in Figure 1b, the modal composition of the granite sample was calculated. This granite sample consists of quartz (30.1%), plagioclase (37.3%), K-feldspar (23.3%), and biotite (9.2%). The modal composition indicates that the granite is a monzogranite [19]. The modal composition was calculated, considering the composition of the abovementioned four minerals, but only without counting the porosity of the granite sample. In addition, the value of porosity (void% in rock) was needed for our simulation work and we obtained that information from the literature, which will be discussed in Section 2.3.

#### *2.2. Experimental*

The electric current–voltage (*I*–*V*) properties of our granite sample were measured. Granite is a well-known dielectric material [3–6] whose resistivity is quite large (1010 <sup>Ω</sup>·m [20]). Thus, when measuring the electric properties of granite, the surface leakage current should be avoided. The sample was placed in a resistivity chamber (Model 24, ADCMT, Saitama, Japan) where the leakage current of the sample was eliminated by a guard electrode of the chamber. The diameter of the upper electrode of this chamber was 6 mm. An ultra-high resistance meter (5451, ADCMT) was used to measure the *I*–*V* properties of the sample. In order to measure the dielectric properties of the sample, we used an LCR meter (ZM2376, NF, Yokohama, Japan) with the same resistivity chamber. The capacitance *C* and dielectric loss *tanδ* were measured with the range of frequency from 1 Hz to 1 MHz with 500 measurement points of frequency.

(**a**) (**b**)

**Figure 1.** Granite sample. (**a**) Optical microscopic image. (**b**) Micro-X-ray fluorescence spectrometer image identifying quartz (red), plagioclase (green), K-feldspar (blue), and biotite (pink).

#### *2.3. Experimental*

The electric current–voltage (*I*–*V*) properties of our granite sample were measured. Granite is a well-known dielectric material [3–6] whose resistivity is quite large (1010 <sup>Ω</sup>·m [20]). Thus, when measuring the electric properties of granite, the surface leakage current should be avoided. The sample was placed in a resistivity chamber (Model 24, ADCMT, Saitama, Japan) where the leakage current of the sample was eliminated by a guard electrode of the chamber. The diameter of the upper electrode of this chamber was 6 mm. An ultra-high resistance meter (5451, ADCMT) was used to measure the *I*–*V* properties of the sample. In order to measure the dielectric properties of the sample, we used an LCR meter (ZM2376, NF, Yokohama, Japan) with the same resistivity chamber. The capacitance *C* and dielectric loss *tanδ* were measured with the range of frequency from 1 Hz to 1 MHz with 500 measurement points of frequency.

Any kind of rock contains voids (pores), though their volumes vary with minerals inside the rocks. Voids are not seen in the modal compositions of our granite sample from μXRF mapping, but we can assume their proportions from the literature [11]. Voids' electrical characteristics can be understood by the gas discharge phenomena [13]. Air is one of the dielectric materials whose dielectric breakdown under a HV application is a large branch of interest in HV engineering [13,19]. In order to understand and simulate the behavior of voids in rocks, we also measured the *I*–*V* properties of the gas discharge by the same experimental procedure mentioned above.

#### *2.4. Simulation Methods*

In this study, the equivalent circuit model was used for simulation. An equivalent circuit model can simulate any kind of rock, however, here, we will discuss granite. The granite was assumed to be a cube of 1 mm edge length, and the model was composed of 1000 of 0.1 mm × 0.1 mm × 0.1 mm smaller cubes (Figure 2a). The model was created by assigning the minerals and voids to those smaller cubes of a 0.1mm edge length and making up the granite. The modal mineral composition of the granite sample discussed in Section 2.1 was used for the simulation. The porosity (i.e., volume % of voids) varies depending on the type and origin of the minerals as well as rocks. For granite, the literature reported that the porosity was from 0.9% to 2.6% [21,22]. In our work, the equivalent circuit model was created by assuming that the porosity in the granite was 2 vol%. The electrical properties of the granite were calculated by simulating the electrical properties of the minerals in the granite, and the circuit parameters were calculated for the equivalent circuit models. A circuit-preparing software program written in C# was developed to create a circuit file in the Netlist format, and then the circuit simulations were performed by using the LTspice circuit simulator.

**Figure 2.** Concept of equivalent circuit model. (**a**) Granite model and (**b**) equivalent circuit model for a small, divided cube of granite.

The possible combinations of mineral distributions can be thousands of patterns. We generated 10 patterns by using random distributions that can decide which minerals or voids will be placed in among the 1000 cubes (0.1 mm edge length). First, void-equivalent circuit elements were placed in the divided cubes with a probability of 2%, and then mineral-equivalent circuit elements for a void were randomly placed in the remaining cubes with the ratio of minerals determined by our μXRF measurement (Figure 1b), based on the fact that the four minerals considered in granite have very similar specific gravities (i.e., quartz (2.65 [11]), plagioclase (2.56 [11]), K-feldspar (2.63 [11]), biotite (2.7–3.4 [23])). Figure 3 shows an example of a mineral distribution using the software of our simulation program. Figure 3a is the 3D model of the minerals and voids (pores) and Figure 3b is the XY plane image of the granite. As mentioned above, the placement of the minerals and voids was selected randomly by our software. When the distribution of the minerals and pores/voids were settled by our software, an equivalent circuit was created automatically to represent the granite sample with the determined distributions of minerals and voids.

**Figure 3.** Example of mineral distribution in the simulated granite model. (**a**) Example of a granite 3D model, (**b**) XY plane of the granite model shown in (**a**).

#### **3. Results and Discussion**

#### *3.1. Electrical Properties*

Following the instructions of the ultra-high resistance meter, we started to take the data of the *I*–*V* properties from the input voltage of 300 V. The applied voltage was then increased from 300 V to 1000 V with an increment of 100 V. The measurements were carried out with an automatic measurement system developed by the LabVIEW software. The data were taken five times under each voltage and their average values and standard deviations were calculated. The relationship of *J*–*E* was calculated from the measured data of the *I*–*V* properties (see Figure 4). Here, *J* represents the electric current density and *E* represents the electric field calculated from the input voltage divided by the distance between the electrodes (i.e., the thickness of the sample, 1.1 mm). *J* was calculated by dividing the electric current *I* with the surface area with the upper electrode, whose diameter was 6 mm. Figure 4 is a double logarithmic graph. The relationship between *J* and *E* can be expressed in the following equation [13]:

$$J = \sigma E \tag{1}$$

where *σ* is conductivity and the slope of the *J*–*E* relationship indicates the conductivity *σ* of the sample. The average value of conductivity was found to be 9.6 × <sup>10</sup>−<sup>11</sup> S/m and its standard deviation was 9.0 × <sup>10</sup>−<sup>12</sup> S/m. Thus, the value of resistivity (1/*σ*) of our sample was 1.0 × 1010 <sup>Ω</sup>·m, and it was similar to the value found in the literature (1010 <sup>Ω</sup>·m) [20].

**Figure 4.** *J*–*E* relationship of granite sample (measured data) with error bars showing the standard deviation values.

#### *3.2. Dielectric Properties of Granite*

The capacitance *C* and dielectric loss *tanδ* were measured with the range of frequency described in Section 2.1 (1 Hz–1 MHz). The measurements were carried out with 500 measurement points over this frequency range. They were measured five times at each frequency and their average and standard deviation values were calculated. The measured data were plotted in double logarithmic graphs (Figure 5). The bars indicate the standard deviation of the values (average data ± standard deviation value). As shown in Figure 5b, the values of *C* decreased with frequency until 1000 Hz but they became almost constant from 1000 to 10,000 Hz. In the latter range, the average value of *C* was 2.45 pF and its standard deviation was 0.253 pF. There were some peak values found in the *C*–*f* properties. These types of peaks in capacitance *C*, as well as the dielectric constant, were found in the dielectric materials in a different range of frequencies, indicating dielectric relaxation [24]. Dielectric relaxation occurs in composite/polymer materials. Dielectric polarization occurs in dielectric materials when a voltage is applied to them. For an AC source of voltage, the re-orientational motions of dipoles and the translational motions of a charged area differ for different types of materials or polymer chains in a composite. This difference, seen at the re-orientational motion, appeared in a resonance form with some frequencies and, hence, the dielectric relaxation occurred in the composite/polymer materials. It is to be noted, however, that the imaginary part of the permittivity shows the dielectric relaxation more clearly [24], but as it is not a matter of discussion for this paper, we have avoided the details here. Peak values were also found in *tanδ*. The average value of *tanδ* from the range of 1000 to 10,000 Hz was 0.773 and the average standard deviation was 0.151. The dielectric constant *ε* was calculated from the Equation (2) [13,25]:

$$\mathbf{C} = \varepsilon \frac{\mathbf{S}}{L} \tag{2}$$

where *C* is the capacitance, *ε* (*ε* = *ε*r*ε*0) is the dielectric constant, *L* is the thickness of the sample (1.1 mm) and *<sup>S</sup>* is the surface area of the electrode (28.3 × <sup>10</sup>−<sup>6</sup> <sup>m</sup>2). *<sup>E</sup>*<sup>r</sup> is the relative dielectric constant/permittivity of granite (i.e., 4–7 at 100 MHz) [26,27] and ε<sup>0</sup> is the dielectric constant of vacuum/air (8.85 × <sup>10</sup>−<sup>12</sup> F/m). From our measurements, the average value of *ε*<sup>r</sup> was found to be 10.8 at the range of frequency from 1000 to 10,000 Hz (10 kHz). The value was not exactly the same as the literature values (i.e., 4–7 at 100 MHz) [26,27]. However, it should be noted that the frequency used in the above reference was 100 MHz, and in general, the capacitance *C* and the relative dielectric constant *ε*<sup>r</sup> decreased with frequency. Using the calculated values of *ε*, the *ε*-*f* properties of our granite sample were plotted in a double logarithmic graph (Figure 6). The resistance *R* is expressed by the following equation, which was used to calculate the general resistance value [25]:

$$R = \frac{1}{\sigma} \frac{L}{S} \tag{3}$$

where *R* is resistance, *σ* is the electrical conductivity, *L* is length (i.e., 1.1 mm), and *S* is the cross-sectional area (i.e., the surface area of the upper electrode, 28.3 × <sup>10</sup>−<sup>6</sup> <sup>m</sup>2). The problem of an electromagnetic wave in dielectric materials (conductivity *σ* = 0) can be understood by solving Maxwell's equations for electromagnetics [25]. We will not discuss the details here, but by solving Maxwell's equations for dielectric materials, the following relationship can be found [25]:

$$
tan\delta = \frac{\sigma}{\omega\varepsilon}\tag{4}$$

where *tanδ* is the dielectric loss, *ε* is the dielectric constant and *ω* represents the angular frequency (i.e., *ω* = 2*πf* where *f* is frequency). Taking the logarithm of both sides of this equation, we obtain Equation (5) [25],

$$
\log \tan \delta \, \, = \log \sigma - \log \omega - \log \varepsilon \,\tag{5}
$$

This equation reflects the behavior of the dielectric materials under the electromagnetic wave, though it does not reflect the dielectric relaxation phenomenon [25,26]. Due to the polarization in a dielectric material (i.e., granite), the change in the dielectric constant/permittivity depends on the frequency *f* of the applied electric field. There were some peaks in the dielectric properties (i.e., *<sup>C</sup>*: 1.68 × <sup>10</sup>−<sup>12</sup> F (393 Hz), 8.44 × <sup>10</sup>−<sup>13</sup> F (21 kHz), 1.59 × <sup>10</sup>−<sup>12</sup> F (435 kHz), 4.37 × <sup>10</sup>−<sup>13</sup> F (761 kHz)—Figure 5a; *tanδ*:0.109 (18.2 kHz), 0.479 (33.4 kHz), 0.0682 (533 kHz)—Figure 5b; <sup>ε</sup>: 6.54 × <sup>10</sup>−<sup>11</sup> F/m (393 Hz), 3.28 × <sup>10</sup>−<sup>11</sup> F/m (21 kHz), 6.59 × <sup>10</sup>−<sup>11</sup> F/m (435 kHz), 1.70 × <sup>10</sup>−<sup>11</sup> F/m (761 kHz)—Figure 6) at different frequencies in the measured data of the sample. From Equations (2), (4), and (5), it is clear that the values of capacitance *C* and dielectric constant/permittivity *ε* have a negative correlation with the angular frequency (*ω* = 2*πf*)/frequency *f.* The slope of log *ε* at the range of frequency 1 to 1000 Hz was found −1.12 (Figure 6). The theoretical value of this slope is −1 (i.e., Equation (5)) and, thus, our measurement results of the dielectric properties of granite show the same trends as the theoretical analysis of dielectric materials. However, after 1000 Hz of frequency, the values of *ε* started to possess peak values, which indicate the dielectric relaxation in the granite, and as Maxwell's equations on electromagnetics cannot

explain the dielectric relaxation of composite materials, our simulation works would not be able to simulate some of the measured data obtained in that frequency range.

**Figure 5.** Measured dielectric properties of granite with error bars showing the standard deviation values. (**a**) *C*–*f* relationship, (**b**) *tanδ*–*f* relationship.

**Figure 6.** *ε*-*f* properties of granite calculated from the measured results of dielectric properties with error bars showing the standard deviation values.

#### *3.3. Electrical Properties of a Void (Pore in Mineral/Rock)*

Two parallel plate electrodes of aluminum were used to measure the *I–V* properties of the air to understand the electric and dielectric behaviors of voids (pores) in granite. The diameter of the electrodes was 4 mm, and the gap was 0.2 mm for the measurement. An ultra-high resistance meter (R8340A, ADVANTEST, Chiyoda ku, Japan) was used to apply the input voltage and measure the electric current. Figure 7 shows the *I*–*V* properties of the air. The current at low-applied voltages was almost negligible, about 150 pA from 5 V to 669 V. On the other hand, at 670 V the current rapidly increased to 0.1 mA, and at 680 V the dielectric breakdown occurred. Thus, the dielectric-breakdown voltage of air was found from this measurement at 680 V/0.2 mm (3.4 × <sup>10</sup><sup>6</sup> V/m). The dielectric-breakdown voltage of air depends on the air pressure, temperature, etc., and in general, it is known

at 3.0 × 106 V/m (3.0 kV for 1 cm of gap length), thus, our result of air breakdown was very similar to other reports [13,20]. From our measurement results, it is clear that the conductivity of a void (pore) is a voltage-dependent parameter.

**Figure 7.** *I–V* properties of the air (i.e., similar to pores or voids in rocks).

#### *3.4. Equivalent Circuit Models for Minerals*

The electrical conductivity *σ* in minerals depends on the temperature. The relationship between the temperature and electrical conductivity *σ* is expressed by the following equation [11]:

$$
\sigma = A \exp\left(\frac{-B}{k\_b T}\right) \tag{6}
$$

where *A* and *B* are constants, *T* is the absolute temperature, and *kb* is Boltzmann's constant (8.618 × <sup>10</sup>−<sup>5</sup> eV/K). The constants *<sup>A</sup>* and *<sup>B</sup>* for each mineral in the granite, and the relative permittivity *ε*<sup>r</sup> used to calculate the capacitance are shown in Table 1. In this study, we used an absolute temperature of 293.15 K, which is considered as room temperature (i.e., 20 ◦C). As shown in Figure 2b, two resistive elements were placed vertically in one cube, and the value of one resistive element was *R*/2. The resistance values *R* of the minerals were calculated using Equation (3). In this study, the capacitance of a mineral was calculated from Equation (2). Since the value of *C* is the value of the entire mineral in the cube of Figure 2b, the value of *C* for the mineral element was 2*C*. The values of *L* and *S* in Equations (2) and (3) were *<sup>L</sup>* = 1.0 × <sup>10</sup>−<sup>4</sup> m and *<sup>S</sup>* = 1.0 × <sup>10</sup>−<sup>8</sup> <sup>m</sup><sup>2</sup> because the minerals were arranged in 0.1 mm cubes, as shown in Figure 2b.

**Table 1.** Parameters used in this simulation work [9,11].


#### *3.5. Equivalent Circuit Model of Void (Pore in Mineral/Rock)*

A void is a pore in a mineral. In this study, we assumed that the void has the same electrical characteristics as the air. Normally, the air does not conduct current when an electric field is applied. However, when a high voltage is applied, a dielectric breakdown can occur, and a very large breakdown-current flows rapidly. The literature indicated that the breakdown voltage of the air is 3 kV for a 1 mm gap of electrodes and the breakdown current varies up to the scale at kA [13]. Therefore, the resistance of a void (pore) can be simulated by a voltage-dependent resistance (VDR), and the capacitor is substituted to create the equivalent circuit model for voids (pores).

In this study, the granite rock was represented by 1000 finite small cubes of 0.1 mm × 0.1 mm × 0.1 mm, and the minerals were supposed to be placed in the small cubes. From the electric and dielectric properties of the minerals (i.e., the conductivity and dielectric constant/permittivity, respectively), the circuit parameters were calculated. When a void was present in a cube with a 0.1 mm edge length within the granite rock, the breakdown voltage was 340 V/0.1 mm (=640 V/0.2 mm, Section 3.3), as we determined from our experiment on the air-discharge phenomenon. The value of the capacitor used in the equivalent circuit model for a void was calculated from the size of the cube and the dielectric constant *<sup>ε</sup>* of air (i.e., 8.85 × <sup>10</sup>−<sup>12</sup> F/m). As Equation (2) is applicable for any dielectric materials that include a void (pore), this equation was used to calculate the void capacitance that was 0.885 × <sup>10</sup>−<sup>15</sup> F and was used in our simulation work. As we described before, the conductivity of air depends on the input voltage, thus, we used a voltage-dependent resistance (VDR) system which can fit with the measurement data of the *I–V* properties of air. The equivalent circuit was simulated with the LTspice circuit simulator. The *I–V* properties of air (void) were simulated for the input voltage ranging from 1 to 350 V. The *I*–*E* characteristics of the simulation data of a void (pore/air) were prepared and plotted in a double logarithmic graph. The equivalent circuit model of the void and the *I*–*E* characteristics of the void are shown in Figure 8a,b, respectively. The electrical properties of a void (see Figure 8) indicates that before the dielectric breakdown of the air, the electric current through a void can be negligible as the value of the electric current is very low (i.e., 10−<sup>12</sup> A at 1 to 350 V) when comparing to that of near the dielectric-breakdown area (i.e., 10−<sup>3</sup> A at 680 V) of a void. We aimed to make an equivalent circuit model for the void, which can represent the *I–V* properties in a way that could fit the dielectric-breakdown region of the void. Figure 8b compared the *I*–*E* properties of the void obtained via measurement and simulation. From the simulation results, the electric current was 20 pA for the input voltage of 10 V (i.e., 100,000 V/m), and it remained less than 10−<sup>9</sup> A even at 150 V (i.e., 1,500,000 V/m). The values of the electric current were 0.1 × <sup>10</sup>−<sup>9</sup> at 50 V (i.e., 500,000 V/m) and 0.21 × <sup>10</sup>−<sup>9</sup> A at 100 V (1,000,000 V/m)). The electric current increased gradually with the increase in input voltage, as we found it was 15.7 × <sup>10</sup>−<sup>9</sup> A at 200 V (i.e., 2,000,000 V/m) and 0.172 × <sup>10</sup>−<sup>6</sup> A at 300 V (i.e., 3,000,000 V/m). Finally, the current reached 1 mA in between the input voltage of 340 V (3,400,000 V/m) and 350 V (i.e., 35,000,000 V/m). The dielectric-breakdown voltage was found in our measurement at 680 V/0.2 mm (i.e., 340 V/0.1 mm), and, thus, we confirmed the reproduction of the *I–V* characteristics of a void (air) in our simulation work. It is to be noted that, while measuring the dielectric breakdown of the void, the system was very unstable to measure, but from the simulation works, the details of the electric current were calculated.

**Figure 8.** Equivalent circuit model and electrical properties of a void (pore) in minerals. (**a**) Equivalent circuit and (**b**) *I*–*E* properties of a void.

#### *3.6. Simulation Results of Granite Sample*

The patterns of minerals and void distributions can be more than thousands, but to evaluate our simulation works, we created 10 distribution patterns of minerals and voids, and then 10 equivalent circuits of granite samples. In order to evaluate the simulation works, a Ltspice circuit simulator was used to simulate the *I–V* properties of the circuits (i.e., granite) under direct-current (DC) input voltages from 100 to 1000 V. They were also simulated with AC voltage to evaluate the dielectric properties of our simulation process with a vast range of frequencies (i.e., 1 Hz to 1000 kHz).

At first, we will discuss the DC properties of our simulation results. The electric current was calculated for each input voltage for each equivalent circuit (i.e., each mineral distribution in granite). From the 10 simulation results of the circuits, the mean values and standard deviation values of the current were calculated for each input voltage. The average *J*–*E* properties were calculated and plotted in a double logarithmic graph, as seen in Figure 9. The bars indicate the standard deviation of the values (average data ± standard deviation value). We also calculated the conductivity *σ* of the granite sample for both the measurement and simulation results (Figure 10a). The average values of *σ* [pS/m] was 53.5 (measurement) and 36.2 (simulation). The standard deviation of *σ* [pS/m] was 38.5 (measurement) and 2.77 (simulation). The literature value of the resistivity (i.e., the reciprocal of conductivity) of granite was in the order of 1012 <sup>Ω</sup>·cm (i.e., 100 pS/m in the value of conductivity) [21] which is the same as our measured and simulation values. Thus, our simulation results are considered reliable. However, a difference in the conductivity properties of granite was found between the measurement results and simulation works. As shown in Figure 11a, it is clear that the simulation values are almost constant regardless of the input voltages, while the measured values increase with the input voltage. This difference might occur from the Schottky effect observed under a HV applied voltage on dielectric materials [13]. For an electric field *E* at temperature *T*, the emission of electrons/ions can be expressed by the following equation:

$$I = SDT^2 \exp\left(\frac{-e}{k\_b T}\right) \exp\left\{\frac{e}{k\_b T} \left(\frac{eE}{4\pi\varepsilon}\right)^{\frac{1}{2}}\right\} \tag{7}$$

where, *<sup>D</sup>* is the Dushman constant (i.e., 1.20 × <sup>10</sup><sup>6</sup> A/m2 <sup>K</sup>2), *<sup>T</sup>* is the absolute temperature in K, *ϕ* is the activation energy (i.e., similar to a work function in metal), *e* is the elementary charge quantity, and the others were described before. The details will not be discussed here, however, there is a linear relationship between the natural logarithm of electric current *I* and electric field *E*1/2 (Equations (8)–(10)).

$$\text{LnI} = K\_0 E^{\frac{1}{2}} + \ln I\_0 \tag{8}$$

$$K\_0 = (\exp/k\_b T)(\exp/4\pi\varepsilon)^{\frac{1}{T}}\tag{9}$$

$$
\ln I\_0 = \ln \left( SDT^2 \right) - \left( e\phi \right) \tag{10}
$$

**Figure 9.** Average of simulation results of *J*–*E* for granite with error bars showing the standard deviation values.

These relationships can be seen in Figure 10b. For the measured data, the determination coefficient was 0.9991, while the same coefficient was 0.9583 for the simulation results. It is worth noting that the Schottky effect was not considered in this simulation work. The electrical properties of the minerals did not consider any voltage-dependent resistance (VDR), but rather, they used constant electrical and dielectric properties to create an equivalent circuit. It explains why the conductivity did not change with the voltage (Figure 10a). As we used voltage-dependent resistance (VDR) for the voids (pores) that account for only 2% in granite, the voltage change had a very limited effect on the conductivity in the simulated results. However, the mean and standard deviation values of the conductivity were very close to the measured values [21]. In fact, it can be understood from the Schottky equation that the resistivity of the rock changes greatly with the temperature, suggesting that it should be considered in future simulations. Even the recognition of voids (i.e., pores in minerals) and considering them into the equivalent circuit accurately may simulate the HV effect on the rocks more concretely. The other paper by our group discussed this relationship (the *I–V* as a function of temperature) under the HV impulse application on granite [19].

Figure 11 shows the mean values of the capacitance *C* and dielectric loss *tanσ* calculated from our simulation works. The error bars indicate the standard deviation of the values (average data ± standard deviation value). In contrast to the measured values (i.e., Figure 3a), the simulation results showed little frequency dependence. The average value of capacitance *C* was 0.0074 pF in between the frequency region of 1 kHz and 10 kHz. The mean value of the standard deviation in this frequency range was 0.000312 pF. These values are about one-tenth of the measured values (i.e., Figure 5a). The *C*–*f* characteristics of the simulation results were different from the measured values, and the values were constant even in the low-frequency range. However, the standard deviation showed a large variation up to 1 kHz, and then it became completely constant. This result suggests that the trend of the dielectric properties could be reproduced, although the average value of *C* could not be completely reproduced. In contrast to the measured values, the simulation results for the *tanδ*–*f* relationship also showed little frequency dependence. The mean value of the loss factor *tanδ* was 0.079, and the mean value of the standard deviation was 0.032. These values were also about one-tenth of the measured values (i.e., Figure 5b). Figure 5 shows that a dielectric relaxation phenomenon occurs around 1 kHz. In the measured data, the peak values for *C* and *tanδ* were found in frequencies other than 1 kHz, but in the

simulation results, it was not observed at frequencies other than 1 kHz. However, we can confirm that the variation becomes larger in the low-frequency band below 100 Hz. This may correlate with the several peaks that occurred in the measured values. From this result, the dielectric phenomenon was considered to be reproduced, although the average value of *tanδ* could not be reproduced completely. It is worth noting that, unlike the measured data, this simulation work did not consider the dielectric relaxation phenomenon in the consideration, thus, the values of the simulated dielectric properties (i.e., *C*, *tanδ*) differed from the measured data. On the other hand, the general trend of the dielectric properties was similar between the simulation and measurement data.

**Figure 10.** Simulation results and analysis (**a**) simulation results of conductivity of granite, (**b**) Schottky effect on granite at room temperature.

**Figure 11.** Simulation results of dielectric properties with error bars showing the standard deviation values. (**a**) *C*–*f* relationship, (**b**) *tanδ*–*f* relationship.

The dielectric constant was calculated from the simulation results of each of the 10 patterns of equivalent circuit simulation, and the average value was calculated. Figure 12 shows the comparison between the calculated results and the measured data. As shown in Figure 12, the dielectric relaxation occurred in the measured dielectric constant values but was hardly observed in the simulation results. However, the measured and simulated dielectric constant values were close to each other in their magnitude. In addition, in the range of frequency between 1 Hz to 1 kHz, the measured values showed a straight line with a negative slope and a constant value after 1 kHz, while the simulation results showed a slight decrease and a constant value after 1 kHz. Thus, it is clear that, considering Maxwell's equation for the dielectric materials (i.e., most of the rocks and granite studied in this work), the dielectric behavior for a certain range of frequency (1–100 Hz for granite in this work) can be reproduced by our simulation method to some extent.

**Figure 12.** Comparison of dielectric constant with frequency in between the measurement and simulation results.

However, it is to be noted that the equivalent circuit model method has some limitations, as recognized by our simulation results of the dielectric properties and dielectric relaxations. Capacitance values, as well as dielectric parameter/permittivity, are strongly time-dependent (i.e., frequency) parameters under a certain applied voltage. However, the capacitance value of a capacitor in the equivalent circuit model is fixed, as the parameter of capacitance cannot be replaced with a time-dependent parameter in circuit simulator software. Dielectric relaxations can be observed in composite materials with different frequencies of AC properties and are also observed in our dielectric measurements of granite. With our current model, the simulation results could not imitate the dielectric relaxations. However, as the parameter *tanδ* is determined with both dielectric and electric properties (i.e., Equations (4) and (5)), a part of the dielectric relaxation properties of *tanδ* was reproduced by our results. These basic disadvantages can be solved by considering several granite samples with a more accurate knowledge of the compositions and their dielectric properties. Further experiments will be needed to improve this method, which will be the next step in our research interests.

#### **4. Conclusions**

In this study, equivalent circuit models of granite were developed by considering the distributions of minerals and voids in granites. In order to confirm the validity of the equivalent circuit model, we measured the electrical and dielectric characteristics of the granite sample and compared them with the simulation results. The results presented in this article are summarized as follows:


These results suggested that the simulation model provided electrical properties similar to those of the granite sample experimentally measured. We prepared equivalent circuit models to simulate granite, but the models can be applied to any rock in order to understand the electrical and dielectric phenomena in rocks. A model considering both conductivity and the dielectric constant will be useful for understanding the electric pulse comminution. Further research works will be carried out to reproduce the dielectric properties of granite with the equivalent circuit model. We are working on a more equivalent circuit method for different ores to develop this method to be more effective in understanding the electric-pulse application behaviors of different minerals, as well as ores.

**Author Contributions:** Conceptualization, M.K. and A.O.; methodology, K.F., M.K., K.K., N.O., M.F. and A.O.; software, K.F., M.K. and K.K.; validation, K.F., M.K., K.K. and A.O.; formal analysis, K.F., M.K. and K.K.; investigation, K.F., M.K. and K.K.; resources, M.K. and A.O.; data curation, K.F., M.K., K.K. and A.O.; writing—original draft preparation, K.F. and M.K.; writing—review and editing, K.F., M.K. and A.O.; visualization, K.F., M.K. and K.K.; supervision, M.K. and A.O.; project administration, M.K. and A.O.; funding acquisition, M.K. and A.O. All authors have read and agreed to the published version of the manuscript.

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

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** We acknowledge Manabu Mizuhira of Bruker Japan for the elemental map analyses of the sample granite at their application laboratory.

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

#### **References**


### *Article* **Non-Destructive Imaging on Synthesised Nanoparticles**

**Kelvin Elphick 1, Akinobu Yamaguchi 2, Akira Otsuki 3,4,5, Neil Lonio Hayagan <sup>3</sup> and Atsufumi Hirohata 1,\***


**\*** Correspondence: atsufumi.hirohata@york.ac.uk

**Abstract:** Our recently developed non-destructive imaging technique was applied for the characterisation of nanoparticles synthesised by X-ray radiolysis and the sol-gel method. The interfacial conditions between the nanoparticles and the substrates were observed by subtracting images taken by scanning electron microscopy at controlled electron acceleration voltages to allow backscattered electrons to be generated predominantly below and above the interfaces. The interfacial adhesion was found to be dependent on the solution pH used for the particle synthesis or particle suspension preparation, proving the change in the particle formation/deposition processes with pH as anticipated and agreed with the prediction based on the Derjaguin–Landau–Verwey–Overbeek (DLVO) theory. We found that our imaging technique was useful for the characterisation of interfaces hidden by nanoparticles to reveal the formation/deposition mechanism and can be extended to the other types of interfaces.

**Keywords:** scanning electron microscopy; backscattered electrons; electron flight simulation; nanoparticles; synthesis

#### **1. Introduction**

Nanoparticles have been synthesised on metallic electrodes and (non-)conductive substrates. Their properties are known to be controlled by their interfacial structures governed by their formation processes. To date, these interfaces have been predominantly imaged by destructive methods, which can achieve nanometric resolution. As reported earlier [1], the highest resolution can be achieved by (scanning) transmission electron microscopy ((S)TEM) and atom probe imaging. These methods have been used commonly for the nano- to atomic-scale analysis of the junction interfaces. However, they require samples to be milled for electron transparency, introducing possible strain and defects during the sample preparation and hindering the direct correlations between the interfacial structures and electromagnetic properties. Electron beam-induced and -absorbed currents (EBIC and EBAC, respectively) has also been used, especially in semiconductor industries, but they are limited to transport properties with the most conductive layer with a submicron resolution.

On the other hand, our recently developed non-destructive imaging method can be performed by controlling the acceleration voltage in scanning electron microscopy (SEM) without modifying a sample and a device [2]. This method achieves a high in-plane resolution of a few nm without any additional requirements of sample preparation for imaging. By controlling the electron acceleration voltages in SEM, the penetration depth

**Citation:** Elphick, K.; Yamaguchi, A.; Otsuki, A.; Hayagan, N.L.; Hirohata, A. Non-Destructive Imaging on Synthesised Nanoparticles. *Materials* **2021**, *14*, 613. https://doi.org/ 10.3390/ma14030613

Academic Editor: Miguel Monge Received: 18 December 2020 Accepted: 25 January 2021 Published: 29 January 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/).

Chiba 277-8581, Japan

of the electron beam can be manipulated. The corresponding generation of secondary electrons (SEs) and backscattered electrons (BSEs) are generated within the electron plume introduced. Since SEs can be surface sensitive via following scattering processes within the specimen, BSEs are detected in this non-destructive imaging method using an energy filter. Recently, we have demonstrated in situ imaging capability under the current-voltage applications, allowing direct comparisons with the defects and the electrical transport properties [3]. Further, the combinations of spectroscopic and scattering/reflective chemical analysis allowed us to evaluate the origins of the defects, which is ideal as a quality assurance for nano-electronic industries. The defect details and the corresponding transport properties can be fed back to the processes of the device fabrication processes, improving the yields [4].

In this study, we applied our method to the characterisation of nanoparticles. We prepared two types of nanoparticles by X-ray radiolysis and the sol-gel method. By imaging these nanoparticles using our method, clear differences in their interfacial structures were found. They revealed the differences in their formation processes during the synthesis or particle suspension preparation, and confirmed the formation/dispersion models predicted depending on the solution pH used for the particle synthesis or particle suspension preparation. Hence, our imaging method can be highly useful for the understanding of the particle synthesis/dispersion processes and can be fed back to the process optimisation of nanoparticle systems.

#### **2. Synthesis of Nanoparticles and Preparation of Nanoparticle Suspensions**

#### *2.1. X-ray Radiolysis*

X-ray radiolysis was used to synthesise nanoparticles using beam line BL8S2 at the Aichi Synchrotron Radiation Center, Aichi Science & Technology Foundation. As detailed in our previous publication [5], a 100-mL aliquot of 0.37 mol/L (M) Cu(COOCH3)2 (FUJI-FILM Wako Pure Chemical Corporation, Osaka, Wako 1st Grade, Japan) was prepared by diluting the stock solution and mixed with methanol with the volume ratio shown in Table 1. In total, 20 μL of these solutions were spread on Si/SiO2, *n*-Si, or Cu substrate, followed by the exposure of 5 min of synchrotron X-ray radiation to synthesise nanoparticles. X-ray irradiation can generate radicals from the radiolysis of liquids and secondary electron generation from substrates dipped in metallic liquid solution. There are some possible routes for the nucleation, ripping, growth, aggregation, and immobilisation of the particles onto the surface of substrate. In particular, near the substrate surface and the interface between the particles and the substrate, the nucleation, growth, and aggregation of these particles can be controlled by the X-ray irradiation significantly. Therefore, this investigation by non-destructive imaging is significantly worthwhile for the understanding of the physical and chemical mechanisms for the synthesis of particles. Simultaneously, this study can also provide the clue to control the synthesis and immobilisation of the particles.


**Table 1.** List of nanoparticles synthesised by X-ray radiolysis.

#### *2.2. Nanoparticle Synthesis by the Sol-Gel Method and Suspension Preparation*

Monodispersed silica nanoparticles were prepared by using the method proposed by [6]. The average particle radius measured by using TEM images was 280 nm and used in the DLVO (Derjaguin–Landau–Verwey–Overbeek) potential calculation (see Section 3). The stock solution containing synthesised silica particles was washed several times to minimise the salt concentration prior to preparation of the desired silica particle suspensions for investigation with the desired chemical environments. Under the different conditions listed in Table 2, silica particle suspensions were prepared in aqueous salt solution (1 × <sup>10</sup>−<sup>2</sup> <sup>M</sup> KNO3, Sigma-Aldrich (St. Louis, MO, USA), and their pH was adjusted using HNO3 or KOH followed by conditioning the suspensions for 30 min. A tiny volume of each sample was pipetted and deposited on a standard SEM aluminium stub that was left in an oven at 50 ◦C for several hours to let the moisture content evaporate and firmly deposit particles on the stub by capillary forces with the residual moisture, followed by metallization of the stub for the sample conductivity.

**Table 2.** List of the silica nanoparticle suspensions studied.


#### **3. DLVO Potential Calculation**

Potential energy calculation between (a) two silica particles or (b) aluminium stub/plate and a silica particle was performed using the DLVO (Derjaguin–Landau–Verwey–Overbeek) theory, which is a well-known theory for describing the material interactions with the summation of the van der Waals potential (*V*A) and electrical double layer potential (*V*R) [7,8]. If the total potential energy (*V*<sup>T</sup> = *V*<sup>A</sup> + *V*R) is high and positive, particles repel each other; otherwise, particles attract each other. This is a straight-forward theory, which can explain particle coagulation/dispersion in many different colloidal systems, e.g., [9–15]. In our previous study that is relevant to the present study, the DLVO theory was also applied to investigate the particle–particle interactions in the system with the small quantity of water present in agglomeration processes [13]. The following paragraphs will introduce and explain the equations used for the potential energy calculation.

Equations used to calculate the potential energies between similar spherical particles [7,8]:

$$V\_{\rm A} = -\frac{Aa}{12H} \tag{1}$$

$$V\_{\rm R} = \frac{64\pi ankT\gamma^2 \exp(-\kappa H)}{\kappa^2} \tag{2}$$

$$n = \text{ N}\_{\text{A}}\text{C}\tag{3}$$

$$\kappa = \left(\frac{8\pi mz^2e^2}{\varepsilon\varepsilon\_0kT}\right)^{\frac{1}{2}}\tag{4}$$

$$\gamma = \frac{\exp\left(\frac{ze\_z^{\tau}}{2kT}\right) - 1}{\exp\left(\frac{ze\_z^{\tau}}{2kT}\right) + 1} \tag{5}$$

where *A* is the Hamaker (J) constant, *a* is the particle radius (nm), *H* is the inter-particle separation distance, *n* is the number concentration of ions (nm<sup>−</sup>3) defined in Equation (3), *N*<sup>A</sup> is the Avogadro's number (6.022 <sup>×</sup> 1023 mol−1), *<sup>C</sup>* is the concentration of ions (mol/nm3), *<sup>k</sup>* is the Boltzmann constant (1.38 × <sup>10</sup>−<sup>23</sup> J/K), *<sup>T</sup>* is the absolute temperature (K), *<sup>γ</sup>* is the reduced surface potential (unitless), *κ* is the Debye–Huckel reciprocal length (nm−1) defined in Equation (4), *ε* is the dielectric constant of the medium, *ε*<sup>0</sup> is the permittivity of free space (C/Vnm), *z* is the ionic valence, *e* is the elementary charge (C), and *ζ* is the zeta potential (V). The zeta potential values of silica particles [16] and aluminium plate [17] were extracted from the literature and used for the present calculation.

Equations used to calculate the potential energy between plate and spherical particle interactions [18] (in our case, the interaction between the aluminium stub and silica particle):

$$V\_{\rm A} = -\frac{A}{6} \left( \frac{a}{H} + \frac{a}{H + 2a} + \ln \left( \frac{a}{H + 2a} \right) \right) \tag{6}$$

$$V\_{\mathbb{R}} = \frac{128\pi amkT\gamma\_{\text{s}}\gamma\_{p}\exp(-\kappa H)}{\kappa^2} \tag{7}$$

where *γ<sup>s</sup>* and *γ<sup>p</sup>* are the reduced surface potential of the sphere and plate (unitless), respectively.

In this article, the calculated total potential energies were normalized by the thermal fluctuation energy (kT). For the dissimilar plate-particle systems, the Hamaker constant *A*<sup>132</sup> was calculated by using the following Equation [19]:

$$A\_{132} = \left(\sqrt{A\_{11}} - \sqrt{A\_{33}}\right)\left(\sqrt{A\_{22}} - \sqrt{A\_{33}}\right) \tag{8}$$

where *A*<sup>11</sup> is the Hamaker constant of particle 1 in vacuum, *A*<sup>22</sup> is the Hamaker constant of particle 2 in vacuum, and *A*<sup>33</sup> is the Hamaker constant of water in vacuum. These values were obtained from the literature [19,20].

#### **4. Non-Destructive Imaging**

As described in Section 1, the acceleration voltage of the electron beam in SEM was precisely controlled to achieve the corresponding penetration into the layer above and below the buried interface to be investigated. The detailed procedures of the nondestructive imaging we recently developed can be found in [3]. An electron flight simulator, CASINO [21], was used to calculate the number of BSEs to be generated in nanoparticles. For the cupric and silica nanoparticles investigated in this study, the simulations show that BSEs can be generated in the vicinity of the nanoparticle–substrate interfaces by introducing an electron beam accelerated at a series of voltages between 18 and 20 keV and between 8.1 and 8.5 keV, respectively. For the latter case for example, as shown in Figure 1a, BSEs are generated predominantly within the nanoparticles at 8.1 keV, while more BSEs are generated from both the nanoparticles and the substrate at 8.5 keV (see Figure 1b). After the lower-acceleration SEM image is subtracted by the higher-acceleration SEM images, buried interfaces can be revealed.

**Figure 1.** Backscattered electrons generated by the electron beam impacted on the nanoparticles at the acceleration voltages of (**a**) 8.1 and (**b**) 8.5 keV. These histograms are simulated by CASINO program [20].

#### **5. Nanoparticles Synthesised by X-ray Radiolysis**

The nanoparticles, #1 and 6, were imaged as shown in Figure 2a–h, respectively. These images were produced after subtracting two SEM images, which have been taken using different acceleration voltages of 18 and 20 keV. These images need to be aligned, which was carried out by adjusting the positions of the nanoparticles within the orange box shown in each image. The colour changes from magenta to green indicate that there are defects or vacancies within the subtracted image. The magenta colour shown at the edge of the particles in these images is due to its spherical shape as there is no intimate contact between the edge of the particles and the substrate. In addition, the bright and dark regions represent the number of BSEs generated to be more and less, respectively.

Figure 2a shows almost white and bright contrast at the nanoparticle–substrate interfaces, indicating that the interfaces are uniformly formed to generate a sufficient number of BSEs. Some arm-shaped regions appeared in magenta and green colours, indicating BSEs are generated above and below the interfaces, respectively. The magenta- and greencoloured features may indicate that the arm regions of the nanoparticles can be detached by voids. This may suggest that these arms can be formed once the main body of the nanoparticle (the middle region) is formed.

Figure 2b shows the nanoparticles synthesised on a *n*-doped Si(001) substrate with a sheet resistance of 1~10 Ω·cm. The size of the nanoparticles is found to be slightly randomised but maintains elongated shapes as seen in Figure 2a. The nanoparticle–substrate interfaces show broad distributions of contrast. This indicates that some nanoparticles with white bright interfaces are formed in the same manner with those synthesised on the Si/SiO2. However, additional nanoparticles may have moved to form clusters, possibly due to the conductivity of the substrate.

(**a**) (**b**)

(**c**) (**d**)

(**e**) (**f**)

**Figure 2.** Processed images of the nanoparticle–substrate interfaces after subtracting two images taken at the acceleration voltages of 18 and 20 keV on the samples grown on (**a**) Si/SiO2, (**b**) *n*-Si(001), (**c**) *p*-Si, (**d**) *n*-Si(111), (**e**) Ni, (**f**) Al, (**g**) 128◦ Y-cut LiNbO3 and (**h**) PTFE substrates.

By replacing the substrate with *p*-doped Si(001) with a sheet resistance of 1~20 Ω·cm, the clustering of the nanoparticles is slightly suppressed by increasing the separation between the nanoparticles as shown in Figure 2c. The shape also becomes square like. The interfaces stay uniform. Their elongation is recovered by synthesising them on *n*-doped Si(111) as shown in Figure 2d, promoting triangular-shaped particles with closer clustering like those on *n*-Si(001).

On the other hand, the nanoparticles synthesised on the metallic Ni substrate show white bright contrast with magenta colour only without any arms as shown in Figure 2e. The shape and size of the nanoparticles are found to be almost cubic with three-fold symmetry as observed for the Si substrates as described above. Similar structures with more elongation are observed for the Al substrate as shown in Figure 2f. Randomly formed nanoparticles are observed for those synthesised on a 128◦ Y-cut LiNbO3 substrate (see Figure 2g). In addition, some distorted particles are observed as immobilised on the LiNbO3 substrate. They may be due to the chemical interactions with Cu(COOCH3)2 and difference in crystallinity between LiNbO3 and cuprates. They become elongated on Al and polytetrafluoroethylene (PTFE) substrates. These results indicate that secondary electrons from the substrates by the X-ray introduction may contribute to the nucleation, growth, and aggregation of nanoparticles. It should be noted that all these samples maintain consistent interfaces.

Since prominent elongated arm-like features were obtained for those synthesised on Si/SiO2 substrates, we further imaged nanoparticles synthesised under a series of pH between 7 and 9. Figure 3a shows almost white bright contrast at the nanoparticle–substrate interfaces for the Y-shaped nanoparticle, confirming that the interface is uniformly formed to generate a sufficient number of BSEs. There are some minor distributions in colour, where some voids exist at the interface. The inverse Y-shaped nanoparticle is found to be formed on the Y-shaped one as the inverse Y-shaped one has darker contrast, indicating the corresponding interface generates less BSEs, i.e., possible presence of voids.

A similar interfacial structure is observed by increasing pH to 9 as shown in Figure 3b. The contrast between the Y-shaped and inverse Y-shaped nanoparticles become weaker, meaning the interface for the latter one reduces the voids, i.e., the formation of a uniform interface. At the same time, each arm becomes more granular than that for the sample synthesised at pH = 8. This means the nanoparticle with the two overlapping Y-shapes is formed by clustering small circular particles to grow along six crystalline facets.

(**a**)

(**b**)

**Figure 3.** Processed images of the nanoparticle–substrate interfaces after subtracting two images 18 and 20 keV on the samples grown on Si/SiO2 at pH = (**a**) 8, (**b**) 9, and (**c**) 7.

(**c**)

For pH = 7, the growth along the six facets becomes weaker to form less elongated arms as seen in Figure 3c. Again, the bottom Y-shaped nanoparticles have stronger adhesion to the substrate than the inverse Y-shaped one. By increasing pH to 11, the nanoparticles become almost like a sphere by attaching only the bottom centre of them. These results suggest that these cuprates are formed from triangular seed crystals, followed by preferred facet growth along the three directions. These cuprates grow in a Y-shape, whose arm length depends on the substrate and pH, which control the mobility of the seed crystals. By rotating 60◦ the Y-shaped cuprates to overlap each other, two of them can form a hexagonal structure [5]. These results shown here indicate that the composition and crystal shape of the synthesised and immobilised cupric nanoparticles are dependent on the conductivity of the substrates and pH of the liquid solutions. The formation of synthesised crystal can be modified when the relative order of the surface energies is altered or when the crystal growth along certain directions is selectively hindered. These results suggest that the selection of surface crystal structures and the electronic states of substrates play a dominant role in controlling the synthesis of nanoparticles.

#### **6. Nanoparticles Synthesised by the Sol-Gel Method and Their Aqueous Suspensions**

Suspensions prepared by the nanoparticles created by the sol-gel method were imaged as shown in Figure 4a–d, respectively, while some additional images were taken for a sample (i.e., Tug 3) as shown in Figure 5. These images were again produced after subtracting two SEM images, which were taken using different acceleration voltages of 8.1 and 8.5 keV. On the other hand, the results of the DLVO potential calculation are shown in Figure 6.

(**a**)

**Figure 4.** *Cont.*

(**b**)

(**c**)

(**d**)

**Figure 4.** Processed images of the nanoparticle–substrate interfaces after subtracting two images taken at the acceleration voltages of 8.5 and 8.1 keV on the samples of Tug (**a**) 2, (**b**) 3, (**c**) 5, and (**d**) 6.

(**b**)

**Figure 5.** Processed images of the nanoparticle–substrate interfaces after subtracting two images taken at the acceleration voltages of 8.5 and 8.1 keV on the Tug 3 sample with different areas (**a**,**b**).

**Figure 6.** DLVO total potential energies between silica–silica particles with the radius of 280 nm as well as between the Al stub plate and silica particles at pH 2 or 10 and 1 <sup>×</sup> <sup>10</sup>−<sup>2</sup> M KNO3 aqueous solution at 25 ◦C, as a function of interparticle separation distance. The unit of the potential, *V*t is the thermal energy, *k*B*T* [J].

In the processed images shown in Figures 4 and 5, in general, the edges of these particles show bright/white colour, indicating that they generate more BSEs as compared with the other interfaces. This suggests that the particles are pinned by these edges due to the capillary force that is often expressed by the neck shape structure between a particle and a plate in the presence of a small amount of water [22–25]. On the other hand, some edges are shown in green or magenta colour, indicating that they are detached from the substrate due to repulsive force that can be explained by the electrostatic interactions between particles and/or particle and the substrate/stub.

In Figure 4a, some nanoparticles contain grey spots (less BSEs), indicating there are some defects or vacancies formed at the nanoparticle, substrate, or their interface. Minor bright edges are observed while other edges in green or magenta are also seen in Tug 2 (Figure 4a), indicating that the nanoparticles are deposited on the substrate with some instability due to electrostatic repulsion among the highly charged silica particles at pH 10 (Tug 2), as the repulsive DLVO potential interaction among them is shown in Figure 6 (pH 10 SiO2 particle–SiO2 particle; pH 10 Al plate–SiO2 particle). It is more noticeable with Tug 2 prepared at pH 10 (Figure 4a) than Tug 5 at pH 2 (Figure 4c), and this agrees with the DLVO potential calculation shown in Figure 6 (pH 10 SiO2 particle–SiO2 particle vs. pH 2 SiO2 particle–SiO2 particle) explaining the repulsive interaction at pH 10 while the attractive interaction at pH 2 between SiO2 particles. Here, it is worth mentioning that the capillary force can be stronger than the attractive DLVO forces [26] (mainly van der Waals force in our case) in order to keep those particles on the substrate/stub while the repulsive DLVO forces (electrostatic force in our case) can influence the stability of particle deposition on the substrate/stub. Similarly, Figure 4c shows almost uniform interfaces at the middle of the nanoparticles but with some edge defects as shown in the green colour. Tug 5 (Figure 4c, pH 2) has a more flat structure in comparison with Tug 2 (Figure 4a, pH 10) that forms multilayer disposition at the same particle concentration of 0.1 vol.%, indicating more stable adhesion between the silica particles and substrate, as the attractive DLVO potential interaction among them at pH 2 shows in Figure 6 (pH 2 SiO2 particle-SiO2 particle).

Figure 4b shows similar interfacial contrast but with darker regions at the edges of the nanoparticles whose suspension was prepared at pH 10 and 0.01 vol.% silica. Furthermore, some interparticle spots in green colour are observed. This indicates that the silica nanoparticles do not have a perfect spherical shape or have repulsive interactions to separate the nanoparticles from each other. On the other hand, Figure 4d shows no interparticle spots in the green colour, indicating attractive interactions between particles as agreed with the DLVO potential calculation (Figure 6, pH 2 SiO2 particle–SiO2 particle). The more uniform colour regions are observed in the Tug 3 (Figure 4b) and Tug 6 (Figure 4d) prepared at 0.01 vol.% silica, where there is a higher amount of moisture content than 0.1 vol.% (Tug 2, Figure 4a and Tug 5, Figure 4c), and it indicated that silica nanoparticles are firmly deposited on the stub by the capillary force with the residual water.

In terms of the interactions between the aluminium stub and silica particles, they are more repulsive at pH 10 and 0.1 vol.% silica particles (Tug 2) than at pH 2 and 0.1 vol.% silica particles (Tug 5), as shown in Figure 4a,c and agreed with the highly positive DLVO potential energies in the former sample condition (Figure 6, pH 10 vs. pH 2) and literature [16]. It can be also explained by the isoelectric point (IEP) of silica, which is around pH 2, where silica particles can coagulate while pH 10 is where the solution pH is far away from IEP and thus silica particles repel each other [16]. Comparing between 0.1 vol.% (Tug 2, pH 10) and 0.01 vol.% (Tug 3, pH 10), the former solid concentration provides more repulsive–instable interactions. It can be explained by the presence of a higher number of particles that are deposited on the aluminium stub with the repulsive nature of interactions.

SEM images were also taken in the other three different areas for each sample to confirm the representativity of the observed images shown in Figure 5 that are the processed images of Tug 3 (pH 10, 0.01 vol.%) as an example. No apparent difference is observed for all the four samples, apart from a minor difference in the particle orientation for the Tug 3. It can be explained by a non-even particle distribution as shown in Figure 5.

#### **7. Summary**

By using our non-destructive imaging method, we imaged nanoparticles synthesised by X-ray radiolysis and the sol-gel method. The X-ray radiolysis is found to initiate the formation of a triangular seed crystal, followed by growth along three facet directions. The sol-gel method, on the other hand, forms spherical nanoparticles, which are pinned to the substrate at the interface and clustered randomly. These crystallisation, deposition, and aggregation processes can be controlled by the substrates, pH, and density as expected and agreed with colloidal DLVO theory. Our imaging method can offer an ideal feedback to achieve precise control of the synthesis processes.

**Author Contributions:** All authors contributed to write this article. A.Y. synthesised nanoparticles and A.O. prepared nanoparticle suspensions. K.E. made all the imaging. A.H. developed the imaging method and analysed the images with K.E., A.O. and N.L.H. performed DLVO potential calculation. All authors have read and agreed to the published version of the manuscript.

**Funding:** The imaging work is partially supported by Japan Society for the Promotion of Science (JSPS)-Engineering and Physical Sciences Research Council (EPSRC) Core-to-core programme (EP/M02458X/1) and Japan Science and Technology Agency (JST) Core Research for Evolutional Science and Technology (CREST) (JPMJCR17J5). The authors thank financial and technical support by JEOL UK to develop the non-destructive imaging. The nanoparticle synthesis and characterisation were partially supported by European Soft Matter Infrastructure (EUSMI) and the visiting professorship of the Institute for Solid State Physics in the University of Tokyo, respectively.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data is contained within the article and available on request with following the guideline set by the University of York (UK).

**Acknowledgments:** A.O. wishes to acknowledge his thanks to Schofield for the silica particle synthesis and Hamane for TEM imaging used for size measurement of silica nanoparticles.

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

#### **References**


*Article*
