**Improving Conductivity in Nano-Conduit Flows by Using Thermal Pulse-Induced Brownian Motion: A Spectral Impulse Intensity Approach**

#### **Ugur Tuzun 1,2**


Received: 4 August 2019; Accepted: 8 September 2019; Published: 17 September 2019

**Featured Application: Nanoparticle suspensions are used in a variety of applications involving the transport of reagents and products to and from structured material surfaces. The more e**ffi**cient alignment of particle layers adjacent to, and in contact with, a reactive surface, such as found in a fuel cell, bio-medical and electrochemical device, is often reliant upon the e**ff**ective control of the dynamics of particle assembly in narrow conduits. One such control is possible by spectral thermal pulsing to generate a controlled Brownian motion of particles. This is demonstrated by numerical simulation here to be highly e**ff**ective when particles are conveyed in viscous fluids close to a neutrally buoyant condition.**

**Abstract:** Inter-particle and particle-wall connectivity in suspension flow has profound effects on thermal and electrical conductivity. The spectral impulse generation and the imparting of kinetic energy on the particles is shown through a mathematical analysis to be effective as a means of achieving an approximate equivalent of a Langevin thermostat. However, with dilute suspensions, the quadratic form of the thermal pulse spectra is modified with a damping coefficient to achieve the desired Langevin value. With the dense suspension system, the relaxation time is calculated from the non-linear differential equation, and the fluid properties were supported by the viscosity coefficient. A "smoothed" pulse is used for each time-step of the flow simulation to take care of the near-neighbor interactions of the adjacent particles. An approximate optimal thermostat is achieved when the number of extra pulses introduced within each time step is found to be nearly equal to the co-ordination number of each particle within the assembly. Furthermore, the ratio of the particle kinetic energy and the thermal energy imparted is found to be never quite equal to unity, as they both depend upon the finite values of the pulse duration and the relaxation time.

**Keywords:** nanoparticle suspension; Brownian motion; spectral thermal pulsing; DEM simulations; Nano-device applications

#### **1. Introduction**

The shallow flows of nano-particles in planar conduits is a common feature of many industrial device applications, as well as being one of the primary generic routes used for transporting reagents and products to and from structured material surfaces as part of selective heterogeneous surface interactions, ranging from adsorptive separations to reaction catalysis; see for example (Pawar and Lee, 2015) [1], (Merino-Garcia et al. 2018) [2].

The more efficient alignment of particle layers adjacent and in contact with a reactive surface, such as found in fuel cell, bio-catalysis and electrochemical surface coating applications, is profoundly important to the efficiency of the processes required to generate, store and release renewable energy, as well as providing uniform surface scaffolds for varied environmental emission detection and separation

devices. Previous work using pore-scale CFD simulations (see for example Crevacore et al. 2017) [3] has investigated the effects of varying the direction of the acceleration due to gravity on colloidal suspensions from orthogonal to parallel to the fluid flow direction, and the resulting impact on Brownian motion and particle settling.

Albert Einstein (1956) [4] formulated the mathematical framework for the kinetics of the Brownian motion of particles in suspension by considering the auto-correlation of fluctuating velocities and the spectral density of the intensity of fluctuations. It is possible to augment the physics of pressure-driven dense assembly flows of near-neutrally buoyant particles by combining the effects of (i) the lubrication forces due to the conveying fluid-immersed particle interactions and (ii) the drag forces due to fluid viscosity with those due to (iii) the thermal pulse-induced Brownian motion of particles.

In a series of DEM simulations of sub-micron Poiseuille flow in a narrow planar channel, the amelioration effect of the quadratic pulse intensity is demonstrated for different combinations of pulse intensity, relaxation time and time-step of simulations as a function of the solid fraction of the suspensions. Most significantly, it is possible to re-direct the clustering of particles in better alignment with the conduit wall surface when thermal pulses are applied. In the absence of pulses, the pressure-driven flows tend to favor a shear-induced migration of particles towards the center of the conduit, where the bulk shear gradient is reduced to zero due to symmetry; see (Davis et al. 2008) [5].

#### **2. Modeling of Lubrication and Viscous Drag Forces**

With dense particle suspensions, the forces between particle pairs are dominated by the lubrication limit, which holds when the gap width between particles, h, is much smaller than their mean diameter, . *h*. The interaction needs to be modified to take account of the roughness of the particle surfaces (no matter how small this may be) to avoid a singularity at h = 0. The interaction needs further augmentation to indicate what happens at h = 0, and for this eventuality a simple collision rheology is introduced as follows.

When two perfectly smooth particles in a fluid in slow flow come close together, the fluid mechanics may be approximated by the lubrication limit. In this limit, the interactive force for approach or departure normal to the smooth surfaces at the surface-to-surface distance h is proportional to . *h*/h. It is therefore never possible for perfectly smooth particles to have a solid contact (h = 0) at a finite velocity. In nature, physical surfaces are never perfectly smooth, and the interaction deviates from the ideal. A roughness dimension has to be introduced. In a paper by Jenkins and Koenders (2005) [6], it is argued that the rough surface effectively behaves as a permeable medium, so that the fluid in squeeze flow can escape by flowing parallel to the surface through the asperities on the surface. In this case a finite relative velocity is possible, and the interactive force behaves more like . *h*/(h + h0), where h0 is a measure of the surface asperity size. Now it is possible to reach h = 0; the particles can touch.

For the purposes of the investigation, suspension flow is here approximated by the motion of an aggregate of neutrally buoyant particles of more or less equal diameter. Then the lubrication force F12 is calculated by using the following equation;

$$\mathbf{F\_{12}} = -3\pi\eta\mathbf{D}^2/8 \left(\mathbf{h} + \mathbf{h\_0}\right). \left(\left(\mathbf{v\_1} - \mathbf{v\_2}\right).\mathbf{\hat{n}\_{12}}\right)\mathbf{\hat{n}\_{12}}\tag{1}$$

where η is the viscosity, *n*ˆ <sup>12</sup> is the unit normal vector pointing from the surface of particle 1 to the surface of particle 2, *D* is the diameter of the particles, v1 and v2 are the velocities of particles 1 and 2. The force mediated by the fluid in this way works purely normal to the surfaces, and a detailed discussion and derivation is in Vahid (2008) [7]. The tangential interaction is of an order less than the leading term in Equation (1) behaving as . *h*ln (h/*D*), see Jeffrey and Onishi (1984) [8].

When particles touch at a finite velocity, a collision mechanism is invoked. The normal velocity is updated using the coefficient of restitution, e using the unit normal vector

$$(\mathbf{v}\_2 - \mathbf{v}\_1)^/. \,\hat{\mathfrak{n}}\_{12} = -(\mathbf{1} - \mathbf{e}) \begin{pmatrix} \mathbf{v}\_2 - \mathbf{v}\_1 \end{pmatrix}. \,\hat{\mathfrak{n}}\_{12} \tag{2}$$

*n*<sup>12</sup> is introduced as the tangential vector such that *n*.*n* = 0. The tangential velocity is unaffected by the collision, so

$$(\mathbf{v}\_2 - \mathbf{v}\_1)^/. \,\overline{n}\_{12} = \begin{pmatrix} \mathbf{v}\_2 \ - \ \mathbf{v}\_1 \end{pmatrix}. \,\overline{n}\_{12} \tag{3}$$

and momentum is conserved,

$$(\mathbf{m}\_1 \mathbf{v}\_1 + \mathbf{m}\_2 \mathbf{v}\_2) = \mathbf{m}\_1 \mathbf{v}\_2 + \mathbf{m}\_2 \mathbf{v}\_2 \tag{4}$$

These equations have been solved (with a symbolic manipulation program) for v1 / and v2 / and then implemented in the computer program.

In order to characterize the viscous drag forces in a particle suspension flow, two alternative approaches may be introduced: The "ensemble-averaged" option, based on a continuum consideration of the interstitial fluid, is to choose the width of the channel and the "suspension viscosity" that is the viscosity of the particle-fluid mixture. In general, the mixture viscosity is a solid fraction dependent factor greater than that of the interstitial fluid, see for example Thomas (1965) [9]. The second option is to choose the particle diameter, and the viscosity of the interstitial fluid thus then allows the calculation of a fluid-drag coefficient as a function of particle size and sphericity; see for example Seville et al. (1997) [10].

For suspension, Reynolds numbers close to Re = 1 (i.e., Stokes dynamics) of the suspension comprising near-neutrally buoyant particles, the effect of the fluid viscosity on the collisional dynamics of the particles can be estimated with the "ensemble-averaged" approach. One such estimate is provided by Davis et al. (2008) [5] where the following formula is derived for the shear viscosity of an isotropic medium;

$$\overline{\eta} = \begin{bmatrix} \mathfrak{A}\eta \phi \mathrm{N}\_{\mathbb{C}} \,/\, 40 \big| \, \mathrm{(D/h)} \, & \,\mathrm{(D/h)} \,\tag{5}$$

Here, η is the suspension viscosity, Nc is the average number of nearby neighbors (Nc = 6 is a good average for a dense suspension) and D/*h* is the ratio of particle diameter to mean surface-to-surface distance, and is a solid fraction-dependent quantity that was estimated by Torquato et al. (1990) [11];

$$\begin{aligned} \text{i}\,\tilde{h}/\text{D} &= (1-\phi)^3/12 \,\phi \,(2-\phi) \end{aligned} \tag{6}$$

where the value of the solid fraction, φ typically ranges from 0.1 to about 0.5 through aggregation of dilute to dense suspensions; see for example Clayton et al. (2012) [12].

#### **3. Characterization of Dem Simulation Results of Pressure-Driven Suspension Flows**

The simulation runs for the external pressure-driven particulate suspension flow case are most easily characterized by means of the Péclet number (see also Morris and Boulay (1999) [13] and Frank, et al. (2003) [14], which is defined as

$$P\_t = \frac{6\pi\eta\dot{\eta}^3}{k\_B T} \,\overline{r}^3 \tag{7}$$

where <sup>η</sup> is the fluid viscosity, . γ is the shear rate, *r* is the mean radius of the particles, *kB* is the Boltzmann constant and *T* is the thermodynamic temperature.

In order to determine these parameters, the graphs of the profiles generated by the simulations are used. For . γ, the velocity profile is employed and an average is calculated from the maximum and minimum velocities and the width of the channel Hp;

$$
\dot{\gamma} \not\cong \mathcal{Q}(\mathbf{v}\_{\max} - \mathbf{v}\_{\min}) / \mathcal{H}\_p \tag{8}
$$

*kBT* is estimated from the equi-partition assumption and is found from the value of the kinetic energy profile at the center of the channel.

A number of 3D DEM simulation results is presented, both with and without Brownian motion for external pressure driven flow depicted by the schematics in section 4 below. It is shown that the simulation captures the migration effect very well, and shows that when the Brownian forces are increased (i.e., reduction of the Péclet number) the migration effect fades. In Figure 1a, the solid fraction profiles are calculated for different values of the Péclet number as a function of the dimensionless distance across the channel yp in terms of the number of particle diameters. When the Brownian motion is switched off; i.e., the Péclet number equals infinity, the case of the near parabolic solid fraction profile is recovered as evidence of the high degree of migration towards the center of the planar channel in Poiseuille flow. In contrast, Figure 1b shows the solid fraction profile for the Péclet number, Pe = 3, evaluated at different time steps for the external pressure driven flow, and it is observed that, as for the very low value of the Péclet number, the solid fraction profile is flatter and the migration effect disappears; see Figure 1b.

**Figure 1.** (**a**) Solid Fraction Profiles across the conduit at different Péclet numbers; (**b)** Solid Fraction Profiles with Péclet number = 3 at different simulation time-steps.

The effects of the presence and absence of the Brownian motion on the solid fraction profiles and the connectivity of the flowing particle assembly by the application of an external pressure for flow show quite clearly that the maintenance of a uniform (flat) solid fraction profile is possible only with a very low Péclet number, which according to Equation (7) also implies that the shear rate and the viscosity are quite low. Hence, in the dilute suspension case depicted here, the increased Brownian motion of the particles tends to encourage particle migration.

These results agree well with the near parabolic solid fraction and particle assembly velocity profiles reported earlier by Morris and Boulay (1999) [13] and Frank et al. (2003) [14]. Statistics of microscopic data, such as the particulate fluctuation energy, can also be presented from DEM simulation results; see Ibrahim (2012) [15].

#### **4. Implementation of Brownian Motion by Thermal Pulses**

For small particles in the sub-micron range, Brownian agitation introduces extra fluctuations, but this motion is thermally induced. While shear-induced fluctuations are position-dependent in flow in a conduit, the thermal fluctuations are position-independent (assuming an isothermal set-up), and therefore, when thermal fluctuations are greater or of the same order as the shear-induced ones, the migration pattern previously observed with shear-induced particle motion is diminished by this effect.

Brownian motion is implemented through a time-dependent, fluctuating force on the particles. Starting from dilute aggregates, a fluctuating force spectrum is derived, which embodies the properties of the fluid, and the momentum transfer from the fluid which can be carried over to the dense suspension case; see Ibrahim (2012) [15] for details. By adding fluctuations, the aggregate increases its

energy, which has to be conducted away through the walls of the flow channel. The thermal conduction process is initiated by giving the walls a fixed temperature.

A one-dimensional problem is examined in the presence of thermal fluctuations, which were introduced to the simulation by the spectral intensity approach. The problem is studied for both cases, dilute and dense slurries. The results for the dilute case were compared with the Langevin method; see for example Davidchack et al. (2009) [16] and Koenders et al. (2012) [17]. Here, the quadratic velocity gives the same answer, while the quadratic displacement is not adequately reproduced. Therefore the quadratic equation is modified in the properties of the fluid in terms of the damping coefficient, and the Langevin value is achieved for small values of the pulse duration compared to the relaxation time; see Section 4.1 below.

The same exercise is done for the dense system where the relaxation is calculated from the non-linear differential equation, and the fluid properties were supported by the viscosity coefficient. It was also observed that the ratio of the kinetic energy and the thermal energy never equal to unity, as they depend on the ratio of the pulse duration and the relaxation time.

The recent work by Ishizuka et al. (2016) [18] also has demonstrated the use of pulsed gas flow in controlling the solids mass flow rate in the riser and downer sections of a triple-bed circulating fluidized reactor. The modulations of pulse width and pulse density (frequency) are used to optimize solids circulation rates by a combination of a low-pass filter and the electrical pulse voltage facilitated by a switching power supply.

#### *4.1. Spectral Intensity of Thermal Energy Fluctuations*

In order to introduce thermal fluctuation to the suspension flow simulation, a spectral intensity approach of impulse distribution has been used on particles for a given fluid viscosity and temperature. The mean quadratic pulse intensity is dependent on the pulse duration and relaxation time, which is calculated by an auto-correlation function, whilst assuming that the pulse duration λ is much smaller than the relaxation time.

The mean quadratic impulse value is given by;

$$<\mathbf{I}^2> = \int \mathbf{I}^2 p(\mathbf{I}) d^2 \mathbf{I} \ = -\frac{\beta}{\pi} \frac{\partial}{\partial \beta} \int e^{\beta \mathbf{I}^2} d^2 \mathbf{I} \ = 1/\beta \tag{9}$$

The value of β is proportional to the mean kinetic energy imparted where β = *m*/ - 2*b*<sup>2</sup> *kB T* identified in two-dimensions for small values of the pulse duration, λ by replacing the mean kinetic energy by *kBT*, instead of 1/2 *kBT*, (which is the one-dimensional value).

It also follows that b = 2 2*m*<sup>γ</sup> <sup>λ</sup> , where *m*/γ is the mean relaxation time. Here, the pulse duration, λ << *m*/γ.

To generate the impulse distribution, the angle is chosen to be entirely random. So, if the probability to encounter a value of I with components between I1 and I1 + *d*I1 and I2 and I2 + *d*I2 is

$$\mathbf{p}(\mathbf{I}\_1, \mathbf{I}\_2) = \frac{\beta}{\pi} e^{-\beta l^2} d\mathbf{I}\_1 d\mathbf{I}\_2. \tag{10a}$$

then

$$\mathbf{p}(\mathbf{I}, \theta) \;= \; \frac{\beta}{\pi} e^{-\beta \mathbf{I}^2} \, \mathbf{I} d\mathbf{I}\_1 d\theta \tag{10b}$$

Integrating over the angle then gives the distribution of the magnitude of the impulses yields;

$$\mathbf{p}(\mathbf{l}) = 2 \frac{\beta}{\pi} e^{-\beta l^2} \mathbf{l} \mathbf{d} \mathbf{l} \tag{11}$$

In practice, the impulses at very high values of I cannot be attained. In order to control this limitation, consider the values generated up to a maximum Imax. Integrating the frequency distribution in the interval 0 < I < Imax results in

$$\int\_{0}^{u} P(\mathbf{I}) \, \, = \int\_{0}^{u} 2\beta \mathbf{I} e^{-\beta \mathbf{I}^{2}} \, \, = \, 1 \, - \, e^{b\mathbf{I}\_{\text{max}}^{2}} \, \tag{12}$$

The value of Imax is chosen such that 99% of all values are covered. In practice values are stored in an array, according to frequencies that correspond to Equation (11), up to the value of Imax, which occurs only once. The angle is chosen at random between 0 and 2π.

The simulation results can then be analyzed in terms of three key parameters: Pulse duration, relaxation time and simulation time step for particle motion.

#### *4.2. Spectral Intensity Generation of Thermal Impulses in 1-D and 2-D Conduit Flows*

Using the framework, which has been developed first for one-dimensional systems, the effects of thermal agitation are also introduced as a fluctuating force in a two-dimensional setting. Thus, the motion of a dense particle slurry that is driven through a conduit by a pressure, combined with Brownian motion, can be studied. The applications for such a system are obviously in nano-technology and in the separation and manipulation of systems in Chemical Engineering where small, sub-micron sized particles are present.

The simplest simulation is an aggregate in a box that is agitated by thermal motion; see Figure 2 below. There are fixed walls at the top and bottom, and periodic boundary conditions left and right. Initially, the wall particles in Figure 2 are given the same fluctuating velocity as the bulk particles; that is, the thermal motion is imparted to them. This allows them to acquire velocity whilst they are kept in the same position. Each particle in the bulk is given an impulse force from the distribution in Equation (11). The resulting particle distribution in a box of 200 particles is shown in Figure 3a.

**Figure 2.** Simulation Channel for Particle-Fluid Suspension.

**Figure 3.** (**a**) Particle Distribution after t = 50,000 time-steps with *kBT* = 3.5 <sup>×</sup>10<sup>−</sup>3; (**b**) Particle Distribution with 'Cold' Walls with *kBT* -> 0 (t > 50,000).

If the walls are left completely "cold", the effect disappears, but a significant migration away from the walls appears as seen in Figure 3b. In order to experiment with the fluctuation velocity which the wall particles attain, the mass of the wall particles is varied. It transpires that for a mass ratio of about 20/1, a uniform (flat) voidage profile is achieved within the conduit. However, this is found to be not sustainable for larger time scales (t > 50,000). The particle migration effect can be suppressed by giving the particles an extra number of pulses in each time-step, so that the effect is a "smoothed" impulse. Further experimentation has established that if the number of extra pulses is equal to the number of nearby interacting particles, an approximate optimum is achieved; see for example; Koenders et al. (2012) [17].

These benchmark simulation experiments are believed to yield two very significant physical outcomes which can facilitate potentially a practical means of controlling the nano-particle agitation and particle separation in conduit flows. The effects suggest that (i) it is possible to impart controlled migration towards the center, and (ii) it is also possible to "clump" particles next to the wall and/or in a periodic cycle to alter the connectivity of the particle bed within the conduit. To achieve these effects in real-life, we would need to impart pulsed impulses to the walls whose materials can be chosen to yield a desired solid density ratio to the particles. A possible application is in the design and operation of electrochemical batteries used for power generation, where the interactions between the electrodes and the charge carrying particle bed can significantly effect the efficiency of electrical charge build-up and controlled release with minimal leakage; see Figure 4 below for a schematic illustration, and refer also to Kim et al. (2013) [19] and Bebelis et al. (2013) [20].

**Figure 4.** Schematic of a Suspension Flow Fuel Cell.

#### **5. Practical Methods to Enhance and Manipulate Brownian Motion in Nano Particle Suspensions**

The literature on nanoparticle suspension conductivity (for example Jang and Choi (2004) [21] and Keblinski et al. (2002) [22]) has considered four possible modes of energy transport in nanoparticle suspension flows: (1) The collisions between carrier fluid molecules, (2) thermal diffusion of nano-particles in suspension, (3) collisions between nanoparticles and (4) thermal interactions between dynamic nanoparticles and the carrier fluid molecules.

The earlier work has found that the temperature-dependent conductivity of nanoparticle suspensions at a specific concentration normalized to the conductivity of the carrier fluid can reveal increases by anything up to two-fold, depending on the fluid temperature. Jang and Choi (2004) [21] present results for Al2O3 particles in water and CuO particles in ethylene glycol, and show that the size of the nanoparticles has a very strong effect on the increase in normalized conductivity achieved. They report that when the fluid viscosity is relatively high like in ethylene glycol, reducing the nanoparticle size to < 10 nm, conductivity increases by two fold; however, when water is used as carrier fluid, the conductivity is reduced with decreasing nanoparticle size.

These results therefore suggest a sensitivity of the ratio of collisions between the nano-particle and carrier fluid molecules and the collisional frequency of the fluid molecules on the values of conductivity achieved. If we also take on board the assertion by Keblinski et al. (2002) [22] that the collisions between nanoparticles are an order of magnitude smaller than the collisions between the nanoparticle and fluid molecules, then the ratio of the pulse intensity to mean relaxation time as a key parameter (see Equation (9)), presented in the thermal pulse intensity spectra analysis presented above, gains further prominence supported by these earlier results reported on the conductivity of nanoparticle suspensions.

Keblinksi et al. (2005) [23] in a more recent study have also shown that the dynamic clustering and aggregation of nanoparticles is more pronounced as the volume concentration of nanoparticles exceeds 5%. This in turn results in a more rapid conductivity due to the increase in particle-particle contacts within the particle clusters. They report results of a several-fold increase in conductivity in consideration of the transformation from randomly closely-packed to loosely-packed spheres; see Keblinksi et al. (2005) [23]. This could only be due to the increased collisional activity of the nanoparticles in loosely-packed particle assemblies providing enhanced connectivity of the nanoparticle assembly.

The pulse intensity spectra analysis presented here shows how the balance between fluid-fluid molecule and fluid-nanoparticle interactions could be altered by affecting the Brownian motion of nanoparticles as function of fluid viscosity and nanoparticle size. Furthermore, it may be possible to make the water-based suspension to behave more like an ethylene glycol suspension, and vice versa, by the imparting of optimal frequency of pulsing for a given concentration and size of nanoparticles. Hence, according to the above analysis, the nanoparticle collisional frequency and its manipulation plays an important role in achieving the optimal value of the ratio of the rate of particle kinetic and thermal energy imparted. This in turn scales the macro-scale conductivity of the nanoparticle suspension relative to the molecular conductivity of the carrier fluid as defined by the well-known Maxwell Model; see Maxwell (1873) [24], which covers the original treatise on electricity and magnetism.

#### **6. Comparison of Large-Scale Simulation Results**

Figures 5 and 6 below compare the thermally pulsed assembly with that of the externally pressurized flow after steady-state is achieved within the flow channel in 60 k simulation time steps with 20,000 particles. The simulations were run for 120 k time steps to achieve accurate and reproducible results. This duration of flow simulation, on average, pumps the whole aggregate through the channel about 40 times, which corresponds to 12 s real-time.

**Figure 5.** Thermally-pulsed assembly flow snapshot at time-step = 100 k.

**Figure 6.** External Pressure driven flow of particulate assembly snapshot at time-step = 100 k.

In both the external-pressure driven flow and thermal-pulsed assembly driven flow scenarios, the suspension is advanced by the creation of periodic cavities adjacent to the planar channel walls coupled with the migration of particles towards the channel center.

There is the significance of the choice of the value of the coefficient of restitution, e = 0.2 (e.g., pharmaceutical compacts; see for example (Bharadwaj et al. 2010) [25] and the fluid viscosity of η = 0.985 Pa-s (e.g., castor oil) with density similar to that of water in the simulation results presented here

to highlight the interactions between highly inelastic "soft" collisions in the presence of moderately high fluid viscosity. These choices were made to allow the assembly simulations to reach "steady-state" within a reasonable computational time frame; see for example, Koenders, et al. (2001) [26].

These results are significant in illustrating the similarity of the nano-conduit flow physics obtained with the creation of Brownian motion in external pressure driven and thermal-pulsed conditions. The advantage of the latter is likely to be realized in applications allowing the energy transfer between the walls and the particle suspension, such as found in fuel cell and electrochemical battery applications. The current DEM simulations, if extended to longer real time flows and larger particle systems, can pave the way to our further understanding of the role of the particle-particle connectivity in particle beds under both Brownian motion-enhanced and -suppressed conditions. The ability to control the extent of Brownian motion in device flows such as seen in Figure 4 above may hold the key to mastering controlled energy storage and delivery with much greater efficiencies achieved to date; see for instance, Mateo et al. (2019) [27].

#### **7. Conclusions**

The results presented of the DEM simulations of neutrally-buoyant sub-micron particles in a flow regime where Re - 0.8 (i.e., Stokes flow), concur well with the earlier observations made of channel flow experiments of suspension flows, whilst highlighting the effects of the controlled Brownian motion of particles via thermal pulsing. The amelioration of the shear-induced migration of particles in planar channel flows is demonstrated by the use of a DEM simulation framework that uses particle collisional dynamics in the presence of the fluid lubrication and viscous drag forces, which also allows for the implicit consideration of the surface roughness of particles through the D/(h + h0) ratio. The simulation results presented in Figures 5 and 6 demonstrate convincingly that the particle-particle and particle-wall thermal and electrical conduction processes are controlled effectively, and could be enhanced to suit practical purposes such as in electrode charging in fuel cells (Canizares et al. 2007) [28] and photovoltaic cells for solar energy storage (Mellor et al. 2016) [29] by manipulating the Brownian motion of the sub-micron particles. Conversely, the Brownian motion could be arrested to a measure to reduce the particle/wall friction experienced in micro-fluidic flows such as in the targeted introduction of bio-pharmaceuticals in the bloodstream, see Pele et al. (2015) [30].

The analytical framework of calculations presented here is capable of incorporation of short-range inter-particle attractive forces and gravitational acceleration to account for particle settling effects, and could be extended to flows in 3-D cylindrical tube configuration; see for example Crevacore et al. 2017) [3]. However, near-planar channel flows (2D) are also of significance in many microfluidic device applications. Most microfluidic devices rely on the operation of a peristaltic pump equivalent, and the thermal pulsing analysis presented here could be adapted to model the effects of pulse width and frequency modulation on collective particle motion in both surface (2D) and tubular (3D) flows; see Ishizuka et al. (2016) [18] for the latter.

**Funding:** Parts of this research received funding from Leverhulme Trust in early stages. Partial latter funding was provided by industrial consultancy activities in the U.S.A.

**Acknowledgments:** This research was supported by a European Federation of Chemical Engineering Energy Section Initiative. Private communications with Professor Manuel Andrés Rodrigo of the Dept. of Chemical Engineering, Faculty of Chemical Sciences and Technologies of the Universidad de Castilla La Mancha, Spain is gratefully acknowledged. Partial external funding was provided by industrial consultancy activities in the U.S.A. The funders had no role in the design of the study, in the collection, analyses and the interpretation of data; in the writing of the manuscript and in the decision to publish the results.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**

1. Pawar, R.C.; Lee, C.S. *Heterogeneous Nano-Composite Photocatalysis for Water Purification*; Micro & Nano Tech. Series; Elsevier Publ.: Amsterdam, The Netherlands, 2015; ISBN 978-0-323-39310-2.


© 2019 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 (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Effect of Soil Reinforcement on Tunnel Deformation as a Result of Stress Relief**

#### **Huasheng Sun \* and Wenbin Sun**

Huaiyin Institute of Technology, 89 North Beijing Road, Huai'an 223001, China; sunwb1969@163.com **\*** Correspondence: sunhuadasheng@126.com; Tel.: +86-18262811276

Received: 23 January 2019; Accepted: 29 March 2019; Published: 4 April 2019

**Abstract:** Adjacent geotechnical engineering activities, such as deep excavation, may adversely affect or even damage adjacent tunnels. Ground reinforcement before excavation may be an effective approach to reduce tunnel heave as a result of stress relief. However, there are few quantitative studies on the effect of soil reinforcement on tunnel deformation. Moreover, the reinforcement mechanism of the reinforced soil and the reinforcement depth are not fully understood. In order to investigate the effect of reinforcing the ground on the tunnel response, a finite element analysis was conducted based on a previously reported centrifugal model test with no ground reinforcement. The effect of the Young's modulus and depth of the reinforced soil on tunnel deformation was analyzed. Soil stresses around the tunnel were also considered to explain the tunnel response. The results revealed that the Young's modulus of the reinforced soil and the reinforcement depth had a significant impact on tunnel deformation as a result of basement excavation. The tunnel heave in the longitudinal direction decreased by 18% and 27% for modulus of the reinforced soil, five times and ten times higher than that of the non-reinforced soil, respectively. The reinforcement depth was effective with regard to controlling the tunnel heave caused by stress relief. This is because the reinforced soil blocked the stress transfer and thus reduced the tunnel heave caused by excavation unloading. It is expected that this study will be useful with regard to taking effective measures and ensuring the safety and serviceability of existing metro tunnels during adjacent excavation.

**Keywords:** basement bottom reinforcement; reinforcement depth; Young's modulus of reinforced soil; tunnel heave; numerical analysis

#### **1. Introduction**

In congested urban areas, the metro plays an important role in the transportation system. The safety and serviceability of existing metro tunnels are always under serious consideration. Owing to the rapid development of underground space in urban areas, the commercial demand for the construction of underground structures in close proximity to metro tunnel lines has been increasing. However, adjacent geotechnical engineering activities, such as deep excavation, may have adverse effects on or even cause damage to nearby existing tunnels. If the induced tunnel deformation and internal forces exceed the capacity of the tunnel structures, segment cracking, leakage, and even longitudinal distortion of the railway track may occur and seriously threaten the smooth travel and safety of the trains in operation. Many studies have investigated the effects of adjacent excavation on existing shield tunnels using various methods, including in situ monitoring [1–3], centrifuge model tests [4–6], numerical analysis [7–15], and semi-analytical methods [16–19]. For example, the main objects of investigation have been the excavation dimension [13,15], relative distance between the tunnel and excavation [4,8,9,11,13], construction and reinforcement methods [10,12], tunnel dimension and physical parameters [9,13], soil density and wall stiffness [6], different constitutive models [14], and influence zone [20].

Chang et al. [2] published an extensive report regarding the damage case history of a shield tunnel as the result of an adjacent deep excavation. In this case history, cracks in segmental linings and the distortion of connected blots were observed. Therefore, it is a major challenge for city designers and geotechnical engineers to evaluate and control shield tunnel responses associated with adjacent excavations. Although various construction and reinforcement methods of controlling tunnel deformation have been proposed by Hu et al. [10] and Liu et al. [12], the existing quantitative research on the evaluation of the bottom reinforcement is insufficient. Moreover, the effect of the soil reinforcement's mechanism on the tunnel deformation has not been clarified.

On the basis of a centrifugal model test reported by Ng et al. [4], the finite element simulation of the studied problem without ground reinforcement has already been done. Ground reinforcement before excavation may be a fresh approach to reduce tunnel heave as a result of stress relief. However, its influence law and deformation mechanism are not clear yet. On the basis of the centrifuge model test without ground reinforcement and the calibrated constitutive models reported by Ng et al. [14], a finite element analysis was conducted to investigate the effect of ground reinforcement on tunnel deformation, both in the longitudinal and in the transverse directions. The effect of the reinforced soil's Young's modulus and the effect of the reinforcement depth under the basement bottom were investigated. The effect of the basement excavation on the deformation of an existing tunnel was analyzed. It is expected that this study may deepen the understanding of tunnel deformation control methods and provide an effective method to control tunnel deformation due to excavation. This study will help in taking effective measures to ensure the safety and serviceability of existing metro tunnels during adjacent excavation for engineering applications.

#### **2. Material and Methods**

#### *2.1. Description of the Simulated Centrifuge Test*

The basis of the numerical simulation was derived from a centrifugal model test conducted by Ng et al. [4]. The tunnel was located just below the center of a pit, as shown in Figures 1 and 2. The centrifuge model test was conducted at the Hong Kong University Science and Technology centrifuge [21–23]. The foundation pit excavation was simulated using the discharge-of-heavy-liquid method, and the excavation depth was controlled using a piezometer. The excavation of the basement was carried out in three steps, and each excavation step was 3 m. Details regarding the process of the centrifuge model tests can be found in the report by Ng et al. [4].

#### *2.2. Finite Element Analysis*

#### 2.2.1. Finite Element Mesh and Boundary Conditions

The finite element program ABAQUS [24] was used to simulate the effect of the basement excavation on the existing tunnel. Figure 3 shows the three-dimensional finite element mesh adopted in the analysis. The mesh dimensions were 1200 mm in length, 990 mm in width, and 750 mm in depth. An eight-node brick element was used to simulate the sand and the diaphragm wall, and a four-node shell element was used to simulate the tunnel. Pin supports were applied to all vertical sides and the base of the mesh to restrain movement in any direction (*x*-, *y*- or *z*-direction).

In all numerical analyses, interface elements were used at the soil–tunnel and soil–basement wall interfaces, unless otherwise specified. Each interface element was a zero-thickness slip element governed by the Coulomb friction law. The friction coefficient (*μ*) and limiting relative displacement (*γ*lim), where the slippage occurred, were controlled by two input parameters for each slip element. The interface friction coefficient *μ* was derived from *μ* = tan *δ*, where *δ* is the interface friction angle, which was taken as 20◦ (i.e., 2/3 of the soil's critical friction angle). The limiting displacement of 5 mm was assumed to achieve the full mobilization of the interface friction.

**Figure 1.** Plan view of the centrifuge model (redraw from Ng et al., 2013) [4].

**Figure 2.** Elevation view of the centrifuge model in model scale (mm) (redraw from Ng et al., 2013) [4].

**Figure 3.** (**a**) Three-dimensional finite element model in model scale (mm) and (**b**) intersection of tunnel and diaphragm wall.

#### 2.2.2. Constitutive Model and Model Parameters

In the experiment, the tunnel and diaphragm wall were made of aluminum alloy and considered elastic materials in the numerical analyses. The Young's modulus and Poisson's ratio of the aluminum alloy were 70 GPa and 0.2, respectively. The sand was simulated using the Mohr–Coulomb model as a first-order model based on the comparison reported by Ng et al. [14]. To describe the mechanical behavior of the soil in this model, five parameters are needed: Young's modulus *E*, Poisson's ratio *ν*, friction angle *ϕ*, cohesion *c*, and dilatancy angle *ψ*. The chosen parameters are listed in Table 1. The parameter selection process and parameter values have been reported by Ng et al. [14].

**Table 1.** Soil parameters used in the Mohr–Coulomb model.


#### 2.2.3. Numerical Modeling Procedure

The numerical modeling procedure is summarized as follows and shown in Figure 4:


**Figure 4.** A sketch of the modeling procedure.

#### 2.2.4. Numerical Modeling Scheme

It should be noted that there is no ground reinforcement in the centrifuge model test. However, the numerical modeling schemes include the cases with and without ground reinforcement. The centrifuge tests provide a benchmark for verifying the correctness and reasonability of the numerical results with no reinforcement. Furthermore, testing of the effect of the reinforcement on tunnel response due to stress relief can be conducted. The ground is assumed to be reinforced by grouting. All parameters of the ground mechanical behavior change after reinforcement may affect tunnel response. A parametric study was conducted to investigate the effect of soil friction angle, cohesion, etc., on the tunnel response. It was found that the strength indexes had only a little influence on tunnel deformation. However, the Young's modulus may play an important role in predicting the tunnel response. Thus, this study focused on the study of the soil's Young's modulus. Figure 5 shows the schematic diagram of the basement bottom reinforcement. The soil mass between the basement bottom and the tunnel was divided into parts 1, 2, and 3, with each part being 10 mm in the model (0.6 m in the prototype). Note that reinforcement S1 was equal to area 1, S2 was equal to area 1 plus area 2, and S3 was the sum of areas 1, 2, and 3. In the first scheme, the effect of the reinforced area's Young's modulus on the behavior of the tunnel was investigated, assuming that the reinforcement area underneath the basement was S1. The reinforced soil's Young's modulus varied between 1, 5, and 10 times the value of the soil's initial Young's modulus. As for the 5-times and 10-times values,

they may not be close to reality but are certainly helpful for practice engineers to judge the effect of soil Young's modulus on the tunnel response. In the second scheme, the Young's modulus of the reinforced soil was constant (*E*<sup>r</sup> = 5*E*s), and the effect of the reinforcement depth on the tunnel's behavior was investigated. The reinforcement depth varied from area S1 = 0.6 m and S2 = 1.2 m to S3 = 1.8 m in the prototype.

**Figure 5.** Reinforcement scheme for basement bottom.

#### **3. Results and Analysis**

#### *3.1. Effect of Reinforced Soil's Young's Modulus on Tunnel Heave*

Figure 6 shows the comparison of the measured and computed normalized tunnel heave and maximum tunnel heave variation with different soil reinforcement moduli, where *H*<sup>e</sup> is the final excavation depth. As can be seen in Figure 6a, the computed values were generally in agreement with the measured data reported by Ng et al. (2013) [4]. The difference between the measurements and the simulation results may be due to the chosen constitutive model. The adopted elastic–plasticity model with Mohr–Coulomb criterion could not capture the strain- and path-dependent soil stiffness. The maximum tunnel heave occurred at the tunnel crown underneath the basement center and decreased gradually as the distance away from the center of the basement increased.

According to the Land Transport Authority of Singapore [25], the maximum tunnel movement should not exceed 15 mm (i.e., 0.17% *H*e). When none of the soil elements were reinforced, the computed maximum tunnel heave was 0.116% *H*e, which overestimated the measured tunnel heave by 50%. In this study, both the computed and measured maximum tunnel heave were within the recommended allowable limit. Because the soil area S1 was reinforced with a Young's modulus that varied between *E*<sup>s</sup> and 5*E*s, the computed tunnel heave decreased greatly from 0.116% *H*<sup>e</sup> to 0.095% *H*e, which overestimated the tunnel heave by 32%. However, the measured tunnel heave was overestimated by 23% as the reinforced soil's Young's modulus increased to 10*E*s. Thus, it is obvious that soil reinforcement underneath the basement reduced the tunnel heave. Further increase of the reinforcement's Young's modulus reduced the tunnel heave, but with a decreasing tendency. For comparison, the overestimated maximum tunnel heave is listed in Table 2. The basic value is the result obtained from the centrifuge model test.

**Figure 6.** Measured and computed (**a**) tunnel heave and (**b**) maximum tunnel heave with different reinforced soil Young's moduli.

**Table 2.** Comparison of computed and measured values.


Note: The −ve sign indicates an underestimation and the +ve sign indicates an overestimation.

#### *3.2. Effect of Reinforced Soil's Young's Modulus on Tunnel Diameter Change*

Figure 7 shows the measured and computed change in the tunnel diameter with different soil reinforcement moduli. The positive and negative values denote the elongation and compression of the tunnel, respectively. It is noted that the described tunnel deformation was located in the tunnel section in a circumferential direction. All three cases predicted that the tunnel lining was vertically elongated and horizontally compressed and that the magnitude of the elongation (Δ*DV*) and compression (Δ*DH*) increased with the excavation depth. All of the computed values underestimated the tunnel diameter change in comparison to the experimental measurements. This can be attributed to the fact that, in the transverse direction, the computed soil stiffness around the tunnel was larger than that in the centrifuge model test. According to the British Tunnelling Society [26], the maximum distortion of a tunnel ((Δ*DV* + Δ*DH*)/*D*) should not exceed 2%. The maximum distortion of the existing tunnel (i.e., 0.09% *D*), which was induced by the basement excavation in this study, was within the abovementioned allowable limit. The increase of the reinforced soil's Young's modulus varied from *E*<sup>s</sup> to 5*E*s, which led to the maximum elongation of the tunnel lining of 9%. Significant changes were not observed when the reinforced soil's Young's modulus increased from 5*E*<sup>s</sup> to 10*E*s. Thus, it was concluded that only a slight influence was exerted by the reinforced soil's Young's modulus on the change of the tunnel diameter.

**Figure 7.** Measured and computed changes in tunnel diameter with different soil reinforcement moduli.

#### *3.3. Effect of Reinforcement Depth on Tunnel Heave*

Figure 8 compares the measured and computed normalized tunnel heave and maximum tunnel heave variation with different reinforcement depths, where *H*e is the final excavation depth. When none of the soil elements were reinforced, the computed maximum tunnel heave was 0.116% *H*e. When the soil area S1 was reinforced, the computed tunnel heave decreased from 0.116% *H*<sup>e</sup> to 0.099% *H*<sup>e</sup> with a constant reinforcement Young's modulus. However, significant changes were not observed when the reinforcement depth increased further. Thus, it is obvious that the soil reinforcement depth underneath the basement reduced the tunnel heave, even for a thin layer. The further increase of the reinforcement depth did not have an effect and was uneconomical.

**Figure 8.** Measured and computed (**a**) tunnel heave and (**b**) maximum tunnel heave with different reinforced soil areas.

#### *3.4. Effect of Reinforcement Depth on Tunnel Diameter Change*

Figure 9 shows the measured and computed changes in the tunnel diameter with different soil reinforcement depths. The positive and negative values denote the elongation and compression of the tunnel, respectively. It is also noted that the described tunnel deformation was located in the tunnel section in a circumferential direction. Tunnel diameter was elongated in the vertical direction and compressed in the horizontal direction after stress relief. The elongated and compressed values increased with the unloading ratio. As the reinforcement depth increased, the elongation and compression of the tunnel lining decreased. Obvious changes were not observed when the reinforcement depth increased further. Thus, it was concluded that the reinforcement depth exerted only a slight influence on the change of the tunnel diameter. When the unloading ratio was 0.5 and the soil area S3 was reinforced, the tunnel diameter change was 0.061%. Compared with the 0.049% change for no reinforcement, a maximum tunnel diameter change ratio of 19% occurred. A summary of the comparison between the computed and measured values is presented in Table 2.

**Figure 9.** Measured and computed change in tunnel diameter with different soil reinforcement areas.

#### *3.5. Analyses of Soil Responses Caused by Excavation*

#### 3.5.1. Stress Distributions of Soil Elements in the Tunnel's Longitudinal Direction

To better understand the tunnel behavior in the transverse and longitudinal directions, the stress distribution of the soil elements around the tunnel are presented in Figures 10 and 11. Figure 10a shows the computed changes in the vertical stress at the tunnel crown in the longitudinal direction for different reinforcement Young's moduli. The positive and negative values denote the increase and decrease in the stress acting on the tunnel lining, respectively. Along the tunnel crown, the vertical stress beneath the basement reduced significantly as a result of the excavation. An almost uniform stress relief was observed beneath the basement. The stress was concentrated beneath the diaphragm wall. This is because the tunnel moved upward due to stress relief, whereas the trend was blocked by the diaphragm wall. The concentrated stress may increase largely beneath the diaphragm wall and then dissipate behind the wall. Thus, it is possible that there will be an increase in stress, which has also been observed in the literature [14]. The vertical stress in the soil increased by more than 100 kPa. After the basement excavation, the maximum change in the vertical stress at the tunnel crown exceeded the allowable limit (i.e., ±20 kPa) set by the Building Department of Hong Kong [27]. In the case where the soil element was not reinforced, the maximum reduction in the vertical stress was 159 kPa. When the reinforced soil's Young's modulus increased from *E*s to 5*E*s, the maximum reduction in the vertical stress was 120 kPa. This was because the reinforced soil blocked the stress transfer and thus reduced the stress release around the tunnel. However, obvious stress changes were not observed when

the reinforced soil's Young's modulus increased from 5*E*s to 10*E*s, which can explain the phenomenon shown in Figure 6. Figure 10b shows the computed changes in the vertical stress at the tunnel crown in the longitudinal direction for different reinforcement depths. As can be seen, the maximum stress reduction was restrained from 159 kPa to 120 kPa when the reinforcement depth changed from 0 to 0.6 m. Further obvious changes were not observed for a reinforcement depth larger than 0.6 m, which is consistent with the tunnel heave change shown in Figure 8.

**Figure 10.** Computed vertical stress change in the soil element at the crown in the circumferential direction for (**a**) different reinforcement Young's moduli and (**b**) different reinforcement depth.

**Figure 11.** Computed earth pressure change in the soil element along the tunnel transverse direction for (**a**) different reinforcement Young's moduli and (**b**) different reinforcement depth.

#### 3.5.2. Stress Distributions of Soil Elements in the Tunnel's Transverse Direction

Figure 11 shows the normal stress change distribution of the soil elements in the transverse direction of the tunnel which was located in the section under the center of the basement. As can be seen, the maximum and minimum stress changes were 80 kPa and almost 0 kPa, respectively, and occurred at the top of the tunnel and at the invert, respectively. It was found that the stress at the invert changed little. This may be due to the rather large stiffness of the tunnel which prevented the stress change of the soil at the tunnel invert. The stress change in the soil elements around the tunnel springlines was 45 kPa. Thus, the tunnel circumference was elongated in the vertical direction and compressed in the horizontal direction, as shown in Figures 7 and 9. As can be seen in Figure 10, in all reinforced cases, only a slight stress change could be observed for the soil elements around the tunnel, in comparison to the non-reinforced case. This can be explained by the diameter change, as shown in Figures 7 and 9.

#### **4. Conclusions**

This study evaluated the effect of the reinforced soil's Young's modulus and the effect of the reinforcement depth on the tunnel response as a result of a basement excavation, based on back-analyzing a centrifuge model test with no reinforcement. The tunnel heave and diameter change and the stress changes both in the longitudinal and in the transverse directions, were analyzed. The obtained results can provide reference for practical engineering in terms of applying construction measures to control tunnel deformation caused by stress relief. Based on the above analyses, the following conclusions were drawn:


**Author Contributions:** Conceptualization, H.S.; methodology, H.S.; software, H.S.; validation, W.S.; formal analysis, H.S. and W.S.; investigation, H.S. and W.S.; resources, W.S.; data curation, H.S.; writing—original draft preparation, H.S.; writing—review and editing, H.S. and W.S.; visualization, H.S.; supervision, H.S.; project administration, H.S.; funding acquisition, H.S.

**Funding:** This study was sponsored by the National Natural Science Foundation of China (grant number 51708245), the Natural Science Foundation of Jiangsu Province (grant number BK20160426), the Construction System Science and Technology Project of Jiangsu Province (grant number 2017ZD124), the Natural Science Foundation for Colleges and Universities in Jiangsu Province (17KJB130003, 17KJA560001), the Project with Science and Technology of Huai'an City (HAG201606).

**Acknowledgments:** We thank Liwen Bianji, Edanz Editing China (www.liwenbianji.cn/ac) for editing the English text of a draft of this manuscript.

**Conflicts of Interest:** The authors declare that there is no conflict of interest regarding the publication of this paper.

**Data Availability Statement:** The computed data used to support the findings of this study are included within the article. Previously reported measured data were used to support this study and are available at [https: //doi.org/10.1139/cgj-2012-0423, https://doi.org/10.1139/cgj-2014-0361]. These prior studies (and datasets) are cited at relevant places within the text as references [4,14].

#### **References**


© 2019 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 (http://creativecommons.org/licenses/by/4.0/).

#### *Article*
