**1. Introduction**

In recent years, there has been significant interest in studying the propagation of elastic waves in solids, as it is crucial for understanding wave propagation mechanisms in engineering applications such as non-destructive testing of structures and the use of new materials [1–8]. Concrete, being a common and popular engineering material, has been extensively studied for its elastic wave propagation mechanism [9–14]. In many concrete structures, circular cavity structures such as pipelines are present, and understanding their response to dynamic loads is important for engineering purposes. At present, there are many different applications for elastic wave propagation in concrete. Many studies have focused on structural damage or defects in concrete. For instance, Ziaja [15] used elastic wave propagation to monitor the state of GFRP-reinforced concrete structural members. They used PZT (lead-zirconate-titanate) sensors to record the changing state of elastic waves caused by cracks and crack propagation in GFRP reinforced concrete structures, considering the material discontinuity caused by cracks and the influence of strain field on wave propagation. Yoon [16] analyzed the applicability of the elastic wave of impact echo (IE) and evaluated six types of prestressed concrete structures using multichannel analysis of surface waves, electromagnetic waves, and shear waves. A more accurate classification method for internal defects in pipelines was proposed by using electromagnetic wave, IE, and principal component analysis (PCA). Beniwal et al. [17] proposed two different ultrasonic imaging techniques designed to use more information contained in the scattered fields for concrete using scattered elastic compression and mode conversion shear wave field modes. Guo [18] established the basic equation of elastic wave propagation in damaged concrete media based on the basic principles of classical elastic dynamics and the damage mechanics model, and derived the fundamental solution of the system. Due to the existence of damage in the structure, the wave response of concrete, including the shape, amplitude, and propagation

**Citation:** Yang, Z.; Sun, C.; Jiang, G.; Song, Y.; Li, X.; Yang, Y. Analytical Study of SH Wave Scattering by a Circular Pipeline in an Inhomogeneous Concrete with Density Variation. *Materials* **2023**, *16*, 3693. https://doi.org/10.3390/ ma16103693

Academic Editor: Yuri Ribakov

Received: 27 March 2023 Revised: 4 May 2023 Accepted: 9 May 2023 Published: 12 May 2023

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

time of ultrasonic waves in the structure, will change obviously. Additionally, many scholars have used the impact-echo method to detect defects in concrete structures [19,20]. Ali [21] described the theoretical basis of a crack detection and location method for concrete samples based on the time for elastic waves generated by crack formation to reach a set of sensors located at the sample boundary, and presented a location method based on acoustic emission detection, and developed a discretization scheme for two-dimensional elastic equations. Uenishi [22] used a two-dimensional in-plane time-harmonic elastodynamics model to analyze the effects of P wave and SV wave incidence on a circular tunnel with lining located at a finite depth from a nearly flat free surface of a homogeneous isotropic linear elastic medium. They also discussed the influence of wavelength and incidence angle, covering layer thickness, and relative compliance on the relative compliance of the linear elastic lining. The results of spalling of lining concrete, buckling of side wall reinforcement, and disengagement of subgrade from invert were given.

In practical engineering structures, the uneven density of concrete can significantly impact the mechanical properties of materials and structures. For example, Lu and Liu [23] analyzed the maximum first principal stress and mid-span deflection increment of density gradient concrete continuous rigid frame bridges under shrinkage and creep effect. Their research aimed to address the problems of excessive mid-span deflection and box girder cracking. The results showed that the effect of shrinkage and creep was reduced by a continuous rigid frame bridge with density gradient concrete. This research provided a theoretical basis for the successful application of continuous rigid-frame bridge with density gradient concrete. Not only do inhomogeneous concrete structures appear in practical projects, but the 3D printing technology of concrete is also becoming increasingly mature, making it possible to prepare concrete materials with gradually functional gradients [24]. Foamed concrete with functional gradients also plays a role in protecting the structure. For example, Strieder [25] used a simplified model to study the influence of gradient concrete material distribution on crack reduction in mass concrete structures. It also demonstrated that graded concrete may help reduce the confinement stress and weaken the risk of cracks during the hardening of concrete. Wang [26] proposed a layered graded foamed concrete-filled tensile honeycomb structure, which achieved multi-level structural protection by adjusting its overall compression deformation mode to layer-bylayer deformation mode. For example, the coagulative density of the new material UR50 ultra-early-strength concrete could reach 2600 kg/m3. This concrete is characterized by ultra-high strength, ultra-high toughness, ultra-impact resistance, and ultra-high durability. Wei et al. [27,28] also researched the mechanical properties and penetration resistance of this high-density concrete. The density of foamed concrete was relatively low, and the density span was large, ranging from 300 kg/m3 to 1600 kg/m3. Therefore, some scholars investigated the performance and structure of foamed concrete with different densities [29]. Foamed concrete backfill has been proven to be an ideal material for improving the bearing capacity of underground engineering structures, except for buried pipelines. Wang [30] studied the effect of foamed concrete backfill on improving the anti-knock performance of buried pipelines. In addition, the study of elastic wave propagation in heterogeneous concrete is of great significance. The propagation characteristics of elastic waves can be used to analyze the structural characteristics of heterogeneous concrete. Metais [31] investigated the impact of multiple scattering on the dispersion curve of phase dry surface waves by considering elastic circular inclusions in an elastic matrix. The dispersion curves were calculated using the global neighborhood algorithm and were inverted to obtain a solution for layered media with linear uniform and isotropic elastic layers. The study quantified the effect of multiple scattering on the results. As the phase velocity of surface waves does not change with frequency, a solution consisting of uniform layers was obtained through inversion. Many scholars have studied the dynamic stress response of homogeneous concrete with defects by numerical simulation or experiment. The dynamic stress response of inhomogeneous concrete with defects is rarely reported. For example, variations in density and shear modulus in inhomogeneous concrete, and structural forms of defects in

concrete. These factors will affect the dynamic stress response of concrete to elastic waves. Even less research has been done on the analytical solutions to such problems. This paper aims to investigate the scattering of SH waves in inhomogeneous concrete containing a circular pipeline using complex function theory and conformal transformation proposed in Refs. [32–36]. In Section 2, the model of density inhomogeneous concrete and the wave field model are introduced. The governing equation is given in Section 3. In Section 4, the stress and displacement fields of the concrete are derived. In Section 5, unknown factors are solved according to boundary conditions, and the expression of dynamic stress concentration factor (DSCF) around the circular pipeline is obtained. In Section 6, the influence of reference wave number and two kinds of density inhomogeneous parameters on DSCF around the circular pipeline is discussed. Finally, Section 7 summarizes the work of this paper.

#### **2. Concrete and Waves Field Model**

#### *2.1. Concrete Model*

The density of concrete will change in the actual engineering structure [23]. In addition, long and thin pipelines such as drainage pipelines will also exist in concrete structures. In view of the above possible situations, this paper assumes that the density of concrete is inhomogeneous. We propose a concrete structure model in which a large volume of concrete contains a relatively small circular pipeline. It is assumed that concrete media is infinite in a two-dimensional plane, and the shear modulus of the concrete *μ* is considered to be a constant value of *μ*0. The density of the concrete is expressed as a polynomial and exponential coupling function, which changes with *x* and *y* in two directions simultaneously and continuously. This problem model is shown in Figure 1. The expressions of concrete density are given by the following equations.

$$
\rho(\mathbf{x}, \mathbf{y}) = \rho\_0 \cdot A(\mathbf{x}, \mathbf{y}) \cdot B(\mathbf{x}), \tag{1}
$$

$$A(\mathbf{x}, \mathbf{y}) = \beta\_1^2 \beta\_2 \, ^2(\mathbf{x}^2 + \mathbf{y}^2),\tag{2}$$

$$B(\mathbf{x}) = \exp(2\beta\_2 \mathbf{x}),\tag{3}$$

where *ρ*0 is the reference density of the concrete, *β*1 and *β*2 are the inhomogeneous parameters of density; The expression of function *A* is the polynomial structure of *x* and *y*, the expression of function *B* is the exponential form.

**Figure 1.** Model of inhomogeneous concrete containing a circular pipeline.

Since the numerical value of the density of the concrete should exist and be real, neither the inhomogeneous parameters *β*1 or *β*2 in the density distribution function expression can be equal to 0. Moreover, the values of inhomogeneous parameters *β*1 and *β*2 affect the variation form of concrete density. The variation of inhomogeneous parameters *β*1 affects the density value in the infinite concrete, while the variation of inhomogeneous parameters *β*2 not only affects the density value, but also directly affects the density distribution in the concrete.

Then, the expression of wave number *k* is given by the following formula:

$$k(\mathbf{x}, y) = k\_0 \sqrt{A(\mathbf{x}, y) \cdot B(\mathbf{x})} \tag{4}$$

$$k\_0 = \frac{\omega\_0}{\sqrt{\frac{\mu\_0}{\rho\_0}}},\tag{5}$$

where *k*0 is the reference wave number, *ω*0 is the circular frequency.

#### *2.2. SH Waves Field Model*

The scattering field model is shown in Figure 1. The radius of the circular pipeline is *R*, and the center of the circular pipeline coincides with the coordinate origin *O*. The density inhomogeneous parameter *β*2 in Figure 1a is a positive number, whereas in Figure 1b, the density inhomogeneous parameter *β*2 is negative. The change from yellow to purple represents that the concrete density changes from large to small. It can be seen from Figure 1 that the density changes in the 2D direction, and is symmetrically distributed along the *x*-axis. Based on the symmetry of density in the concrete, the incident SH wave is assumed to incident horizontally along the *x*-axis. When *β*2 is positive, SH waves are incident from the low density to the high density of the concrete; when *β*2 is negative, SH waves are incident from the high density to the low density of the concrete. The incident direction of the two cases is completely opposite.

#### **3. Governing Equation**

In the Cartesian coordinate system, the wave equation in the concrete with inhomogeneous density is given by the following equation:

$$\frac{\partial^2 \varrho(\mathbf{x}, y)}{\partial \mathbf{x}^2} + \frac{\partial^2 \varrho(\mathbf{x}, y)}{\partial y^2} + k^2(\mathbf{x}, y) \cdot \varrho(\mathbf{x}, y) = 0,\tag{6}$$

where *ϕ*(*x*, *y*) is the displacement in the wave field, which is the function of *x* and *y*.

Based on the complex function theory, a set of complex variables, *z* = *x* + *iy* and *z* = *x* − *iy*, where introduced to transform the wave Equation (6) into the following equation:

$$\frac{\partial^2 \varrho}{\partial z \partial \overline{z}} + \frac{1}{4}k^2(z, \overline{z}) \cdot \varrho = 0,\tag{7}$$

where *k*(*z*, *z*) is expressed as *k*(*z*, *z*) = *k*0*β*1*β*2|*z*| exp[0.5*β*2(*z* + *z*)] in the coordinates of the complex variables.

To solve wave Equation (5), we need to introduce a new set of variables *ζ* and *ζ*

$$\mathcal{L} = w(z) = \beta\_1 (z - \frac{1}{\beta\_2}) \exp(\beta\_2 z) \; , \; \overline{\mathbb{S}} = w(\overline{z}) = \beta\_1 (\overline{z} - \frac{1}{\beta\_2}) \exp(\beta\_2 \overline{z}) . \tag{8}$$

By introducing a new set of variables, the Helmholtz equation with variable coefficients can be transformed into the standard one, allowing for the easy derivation of analytic solutions for displacement and stress fields using the separation of variables method. The standard form of the Helmholtz equation is expressed as:

$$\frac{\partial^2 \varrho}{\partial \zeta \partial \overline{\zeta}} + \frac{1}{4} k\_0^2 \varrho = 0. \tag{9}$$

#### **4. Fields of Displacements and Stresses**

The propagation direction of the incident wave is horizontal and the incident angle is 0◦. In the *ζ*-plane, the displacement of the wave field can be obtained by using the Helmholtz equation in the standard form, and the expression of the incident waves *ϕ*(*i*) is as follows.

$$\!\!\!\!\!\!\!\!\!\/^{(i)}\left(\zeta,\overline{\zeta}\right) = \!\!\!\!\!\!\!\/^{(i}\!\!\!\!\!\/ \_{0}\left(\zeta+\overline{\zeta}\right)\!\!\!\/) \!\!\!\/^{(i)} \tag{10}$$

where *ϕ*0 is the amplitude of the incident wave.

In addition, a circular pipeline exists in an infinite inhomogeneous concrete, and the scattering waves *ϕ*(*s*) from the circular pipeline is,

$$\mathfrak{q}^{(s)}\left(\mathbb{\zeta},\overline{\mathbb{\zeta}}\right) = \sum\_{n=-\infty}^{\infty} \mathbb{C}\_{n} H\_{n}^{(1)}\left(k\_{0}|\mathbb{\zeta}|\right) \left(\frac{\mathbb{\zeta}}{|\mathbb{\zeta}|}\right)^{n},\tag{11}$$

where *Cn* are undetermined coefficients and *H*(1) *n* is the first kind Hankel function of *n*th order.

In the concrete, the displacement field should be the superposition of incident waves and scattering waves displacement, so the displacement fields in the concrete with inhomogeneous density can be expressed as:

$$
\varphi\_t = \varphi^{(i)} + \varphi^{(s)},\tag{12}
$$

where *ϕt* represents the total displacement field.

In the complex plane, the expression of hoop stresses and radial stresses in the concrete with inhomogeneous density is given by the following equation,

$$
\pi\_{rz} = \mu\_0 \left( \frac{\partial \varphi}{\partial z} e^{i\theta} + \frac{\partial \varphi}{\partial \overline{z}} e^{-i\theta} \right) \tag{13}
$$

$$\text{tr}\_{\theta z} = i\mu\_0 \left( \frac{\partial \,\varphi}{\partial z} e^{i\theta} - \frac{\partial \,\varphi}{\partial \overline{z}} e^{-i\theta} \right) \,\text{\textquotedbl{}}\tag{14}$$

where *r*, *θ*, and *z* are the cylindrical coordinates.

By introducing variables *ζ* and *ζ*, Equations (13) and (14) can be transformed into the following,

$$
\pi\_{rz} = \mu\_0 \left( \frac{\partial \varrho}{\partial \zeta} \frac{d\overline{\zeta}}{dz} e^{i\theta} + \frac{\partial \varrho}{\partial \overline{\zeta}} \frac{d\overline{\zeta}}{d\overline{z}} e^{-i\theta} \right),
\tag{15}
$$

$$\pi\_{\theta z} = i\mu\_0 \left( \frac{\partial \,\rho}{\partial \tilde{\zeta}} \frac{d\tilde{\zeta}}{dz} e^{i\theta} - \frac{\partial \,\rho}{\partial \tilde{\zeta}} \frac{d\overline{\tilde{\zeta}}}{d\overline{z}} e^{-i\theta} \right). \tag{16}$$

By substituting the displacement of incident waves and scattering waves into Equations (15) and (16), the stresses field in the infinite inhomogeneous concrete in the *ζ* plane can be obtained,

$$\tau\_{rz}^{(i)} = \frac{1}{2} i \mu\_0 k\_0 \rho\_0 \beta\_1 \beta\_2 \left[ z \exp(\beta\_2 z + i\theta) + \overline{z} \exp(\beta\_2 \overline{z} - i\theta) \right] \exp\left[ \frac{i k\_0}{2} (\zeta + \overline{\zeta}) \right],\tag{17}$$

$$\pi\_{\theta z}^{(i)} = -\frac{1}{2} \mu\_0 k\_0 \rho\_0 \rho\_1 \beta\_2 [z \exp(\beta\_2 z + i\theta) - \overline{z} \exp(\beta\_2 \overline{z} - i\theta)] \exp\left[\frac{ik\_0}{2} \left(\zeta + \overline{\zeta}\right)\right],\tag{18}$$

$$\mathbf{T}\_{\tau\_{\Xi}}^{(s)} = \frac{\mu\_0 k\_0 \theta\_1 \theta\_2}{2} \sum\_{n=-\infty}^{\infty} \mathbb{C}\_{\pi} \left\{ H\_{n-1}^{(1)}(k\_0 |\underline{\zeta}|) \left(\frac{\underline{\zeta}}{|\underline{\zeta}|}\right)^{n-1} z \exp(\beta\_2 z + i\theta) - H\_{n+1}^{(1)}(k\_0 |\underline{\zeta}|) \left(\frac{\underline{\zeta}}{|\underline{\zeta}|}\right)^{n+1} \mathbb{E} \exp(\beta\_2 \overline{z} - i\theta) \right\},\tag{19}$$

$$\pi\_{\theta z}^{(s)} = \frac{i\mu\_0 k\_0 \theta\_1 \theta\_2}{2} \sum\_{n=-\infty}^{\infty} \mathbb{C}\_n \left\{ H\_{n-1}^{(1)}(k\_0|\zeta|) \left(\frac{\zeta}{|\zeta|}\right)^{n-1} z \exp(\beta\_2 z + i\theta) + H\_{n+1}^{(1)}(k\_0|\zeta|) \left(\frac{\zeta}{|\zeta|}\right)^{n+1} \overline{z} \exp(\beta\_2 \overline{z} - i\theta) \right\}. \tag{20}$$

#### **5. Boundary Conditions and Dynamic Stress Concentration Factor of Circular Pipeline**

According to the relationship between the infinite concrete and the circular pipeline, it can be determined that the stress freedom should be satisfied on the boundary of the circular pipeline. Therefore, the radial stresses should be zero on this boundary (|*z*| = *R*) as

$$
\tau\_{rz}^{(i)} = \tau\_{rz}^{(i)} + \tau\_{rz}^{(s)} = 0,\\
$$

According to Equation (21), the radial stresses expressions of the incident and scattered waves are substituted into Equation (19), and the following expressions can be obtained,

$$\sum\_{n=-\infty}^{\infty} \mathbb{C}\_n \mathfrak{f}\_n = \mathfrak{f}\_{\prime} \tag{22}$$

where

$$\xi\_{\boldsymbol{\theta}} = H\_{n-1}^{(1)}(k\_0|\boldsymbol{\zeta}|) \left(\frac{\boldsymbol{\zeta}}{|\boldsymbol{\zeta}|}\right)^{n-1} \boldsymbol{z} \exp(\beta\_2 \boldsymbol{z} + i\boldsymbol{\theta}) - H\_{n+1}^{(1)}(k\_0|\boldsymbol{\zeta}|) \left(\frac{\boldsymbol{\zeta}}{|\boldsymbol{\zeta}|}\right)^{n+1} \boldsymbol{\overline{z}} \exp(\beta\_2 \boldsymbol{\overline{z}} - i\boldsymbol{\theta}), \tag{23}$$

$$\mathbb{E} i k\_0, \qquad \boldsymbol{\underline{\boldsymbol{\zeta}}}$$

$$\mathcal{L} = -i\rho\_0[z\exp(\beta\_2 z + i\theta) + \overline{z}\exp(\beta\_2 \overline{z} - i\theta)]\exp\left[\frac{ik\_0}{2}(\zeta + \overline{\zeta})\right].\tag{24}$$

To solve the unknown coefficient *Cn*, multiply both sides of Equation (22) by *e*−*imθ*e, and perform integration over the interval of 2*π* to obtain Equation (25).

$$\sum\_{n=-\infty}^{\infty} \frac{\mathbb{C}\_{\mathfrak{n}}}{2\pi} \int\_{-\pi}^{\pi} \mathbb{S}\_{\mathfrak{n}} e^{-im\theta} \mathrm{d}\theta = \frac{1}{2\pi} \int\_{-\pi}^{\pi} \mathbb{S} e^{-im\theta} \mathrm{d}\theta,\ (m = 0, \pm 1, \pm 2, \dots, \pm n). \tag{25}$$

Then the dynamic stress concentration factor around the circular pipeline in the concrete can be obtained, which is given by the following formula,

$$
\pi\_{\theta z}^\* = \left| \frac{\pi\_{\theta z}^{(t)}}{\pi\_0} \right|,
\tag{26}
$$

$$\begin{split} \tau\_{\theta\varepsilon}^{(t)} &= -[z \exp(\beta\_{2}\varepsilon + i\theta) - \overline{z} \exp(\beta\_{2}\overline{z} - i\theta)] \exp\left[\frac{i\delta\_{0}}{2}(\xi + \overline{\xi})\right] \\ &+ \frac{i}{\theta\_{0}} \sum\_{n=-\infty}^{\infty} \mathcal{C}\_{n} \left\{ H\_{n-1}^{(1)}(k\_{0}|\xi|) \left(\frac{\overline{\xi}}{|\overline{\xi}|}\right)^{n-1} z \exp(\beta\_{2}\varepsilon + i\theta) + H\_{n+1}^{(1)}(k\_{0}|\xi|) \left(\frac{\overline{\xi}}{|\overline{\xi}|}\right)^{n+1} \overline{z} \exp(\beta\_{2}\overline{z} - i\theta) \right\}^{\prime} \end{split} \tag{27}$$

$$\text{where } \tau\_{\theta\varepsilon}^{(t)} = \tau\_{\theta\varepsilon}^{(i)} + \tau\_{\theta\varepsilon}^{(s)} \cdot \tau\_{0} = \mu\_{0}k\_{0} \varrho\_{0} \beta\_{1} \beta\_{2} \text{.}$$

#### **6. Numerical Results and Discussion**

After establishing the scattering model of elastic waves by circular pipeline in the infinite inhomogeneous concrete, the DSCFs around the circular pipeline can be obtained when SH waves propagate in the horizontal direction. The distribution and variation rule of the DSCFs around the circular pipeline are analyzed and discussed. The dimensionless variables used in the analysis are reference wave number *k*0*R*, density inhomogeneous parameters *β*1 and *β*2, respectively. When other variables are the same as each other and the values of *β*2 are opposite number to each other, the density distribution in the concrete is symmetric about the y-axis. When *β*2 is positive, it can be considered that SH waves are incident horizontally from low density to high density in the concrete. When *β*2 is negative, SH waves are incident horizontally from high density to low density in the concrete.

To analyze the impact of different values of inhomogeneous parameters on DSCF, the values of inhomogeneous parameters adopted in this paper were based on the selected variables in the study of reference [37]. The values of variables in reference [37] can be found in Appendix A. The distribution of DSCFs around the circular pipeline is given in Figure 2 when the density inhomogeneous parameter *β*1 = 0.5. The values of *k*0*R* are 0.5, 1.0 and

2.0, respectively. The density inhomogeneous parameter *β*2 in Figure 2a,b is 0.1 and 0.5. As can be seen from the figure, since the distribution form of the inhomogeneous concrete is symmetric about the *x*-axis, and SH waves are also positively incident along the *x*-axis, the distribution of the DSCFs around the circular pipeline is also symmetric. When the reference wave numbers *k*0*R* become two times larger, the amplitude of the DSCF around the circular pipeline increases. It can be inferred that higher DSCF values will be obtained in the case of high-frequency waves at this density. When *β*2 increases, the maximum value of DSCF also increases, and the extreme point of DSCF shifts toward the back wave surface, and the distribution of DSCF around the circular pipeline becomes more regular. Figure 2c,d shows that when *β*2 is negative, the direction of density distribution in the concrete changes concerning the symmetry of the y-axis. Contrary to the results in Figure 2a,b cases, the amplitude of DSCF decreases with the increase of reference wave number *k*0*R*, and higher DSCF values are obtained in the case of low-frequency waves. Compared with the case when *β*2 is positive, the maximum value of DSCF when *β*2 is negative is smaller. However, in the same way, that *β*2 is positive, the extreme points of DSCF will shift toward the back wave surface of the circular pipeline with the increase of the absolute values of *β*2. Thus, it can be inferred that when SH waves are incident into the concrete from two opposite directions, the distribution law and magnitude of DSCF will be changed.

**Figure 2.** DSCF distribution around the circular pipeline with different reference wave numbers *k*0*R*(*β*1 = 0.5).

By comparing the results of the DSCFs around the circular pipeline in the radial and linear terms of the density given by Jiang [37], the distribution of DSCFs around the circular pipeline given in Figure 2 shows that when the density of the concrete is distributed in the form of a polynomial-exponential coupling, more extreme points appear in the DSCF

curves. The amplitudes of DSCF all appear on the surface, and there is no obvious offset with the increase of reference wave number or inhomogeneous parameters. It shows that when the density of the concrete is distributed in this form, the influence on the distribution of DSCF around the circular pipeline is more severe.

The values of DSCFs under different reference wave numbers *k*0*R* are given in Figures 3 and 4. The values of *k*0*R* are set as 0.5, 1.0 and 2.0, respectively.

**Figure 3.** DSCF distribution around circular pipeline with different reference wave numbers *k*0*R* (*β*2 = 0.5).

**Figure 4.** DSCF distribution around circular pipeline with different reference wave numbers *k*0*R* (*β*2 = −0.5).

Figure 3 shows the distribution of DSCF around the circular pipeline when the *β*2 = 0.5 and *β*1 = 0.3, 0.6, 0.9. The DSCFs around the circular pipeline increase with larger reference wave number *k*0*R*. When *β*1 increases, it can be seen that the distribution of DSCF around the circular pipeline (20~40◦ and 320~340◦) becomes regular, and the number of extreme points of DSCF around here decreases. When *k*0*R* and *β*1 change, the distribution of extreme points of DSCF is dominated by the back wave surface. Figure 4 shows the DSCF distribution around the circular pipeline when the *β*2 = −0.5. It can be seen that when *β*2 is negative, the maximum value of DSCF decreases with the increases of *k*0*R*, but its maximum value decreases slightly. The DSCF around the circular pipeline is complex, and there is no obvious distribution change with the increase of *β*1, and the maximum value of DSCF is much smaller than the result in Figure 3. Thus, it can be found that when the incident direction of SH waves is from low density to high density, the distribution of DSCF in the process is more regular.

In the discussion of the influence of reference wave number *k*0*R* on DSCF around the circular pipeline, it is obvious that the inhomogeneous parameters *β*1 and *β*2 will have a significant influence on DSCF around the circular pipeline. Figure 5 shows the variation of DSCF distribution around the circular pipeline with the values of *β*2, when the *k*0*R* = 0.5, 1.0 and 1.5. We set the *β*1 = 2.0 and changed the *β*2 = 0.2, 0.3 and 0.4. As can be seen from Figure 5, when the *β*2 increases, the maximum value of the DSCF around the circular

pipeline increases along with it. Although the values of *β*2 are small and the increment of each change is small, the amplitude of the DSCFs around the circular pipeline is very obvious, because the change of *β*2 has a drastic impact on the density distribution and value in the concrete. In addition, as the reference wave number *k*0*R* increases into an arithmetic sequence, it is found that the distribution of DSCF near 40◦ and 320◦ of the circular pipeline becomes more regular, and the number of extreme points of DSCF decreases. Figure 6 shows DSCF around the circular pipeline with the values of *β*2 = 0.8, 1.0 and 1.2. When the value of *β*2 is large, the maximum value of the DSCF around the circular pipeline will not change significantly with the values of *k*0*R* variation. Therefore, it can be inferred that when *β*2 increases to a certain value range, the change of reference wave number has a small impact on the amplitude of DSCF around the circular pipeline. In addition with the increase of *β*2, the distribution of DSCF near 20◦ and 340◦ of the circular pipeline gradually becomes complex, and the number of extreme points of DSCF tends to increase. These changes are completely contrary to the changes when *β*2 is small.

**Figure 5.** Distribution of DSCF around circular pipeline with different inhomogeneous parameter *β*2(*β*2 = 0.2, 0.3, 0.4), (*β*1 = 2.0).

**Figure 6.** Distribution of DSCF around circular pipeline with different inhomogeneous parameter *β*2(*β*2 =0.8, 1.0, 1.2), (*β*1 = 2.0).

The distribution of DSCFs around the circular pipeline with *β*1 is given in Figures 7 and 8. The reference wave number *k*0*R* = 0.5, 1.0 and 1.5. Figure 7 set the *β*2 = 0.5, and change the *β*1 = 1.2, 1.6 and 2.0. As can be seen from Figure 7, when the *β*1 increases, the maximum value of DSCF around the circular pipeline increases along with it. When *k*0*R* increases, it can be found that the distribution of DSCF near 30◦ and 330◦ of the circular pipeline changes gradually, and the number of extreme points of DSCF decreases. In addition, *k*0*R* increases in an arithmetic sequence, and the maximum value of DSCF increases in an arithmetic sequence. Figure 8 shows the distribution of DSCF around the circular pipeline at *β*2 = −0.5. SH waves are incident horizontally from high density to low density in the concrete. It can be found that the distribution of DSCF around the circular pipeline is more complex than in Figure 7, in which multiple extreme points appear on the back wave surface, and when *β*1 increases, it

has little influence on the amplitude of DSCF. When *k*0*R* increases, the peak values of DSCF almost do not change, remaining within 0.4, which is smaller than the maximum value of DSCF in Figure 7. Only the distribution of DSCF near 30◦ and 330◦ of the circular pipeline changed slightly. It can be inferred that when SH waves are incident from high density to low density, the changes of *k*0*R* and *β*1 have a very limited impact on the amplitude and distribution of DSCF around the circular pipeline, and no obvious changes will occur.

**Figure 7.** Distribution of DSCF around circular pipeline with different inhomogeneous parameters *β*1 (*β*2 = 0.5).

**Figure 8.** Distribution of DSCF around circular pipeline with different inhomogeneous parameters *β*1 (*β*2 = −0.5).

The distribution of DSCF around a circular pipeline is shown in Figures 2–8, and in most cases, there are significant differences in the distribution of DSCF between 20~40◦ and 320~340◦. In order to analyze and summarize the distribution characteristics of DSCF around the circular pipeline more intuitively, this paper selects three observation points on the circular pipe at *θ* = 20◦, 30◦, and 40◦, and provides the changes in DSCF values under different variable influences in Tables 1–3.

 The results of the values of DSCF at positions *θ* = 20◦, 30◦, and 40◦ on a circular pipeline with *β*1 = 2.0 and *k*0*R* = 0.5, 1.0, 1.5 are shown in Tables 1–3. It can be visually observed that the values of DSCF at the position of *θ* = 30◦ on the circular pipeline are relatively more sensitive to changes in *β*2 and *k*0*R* compared to the other two observation points, based on both horizontal and vertical comparisons of the data in the three tables. It can be inferred that changes in the parameters of inhomogeneous density and reference wavenumber have the most significant impact on the values of DSCF at the position of *θ* = 30◦ on the circular pipeline. Therefore, the position *θ* = 30◦ is selected as the observation point. The continuous variation of DSCF with *β*1, *β*2 and *k*0*R* at the observation point is analyzed and discussed.

**Table 1.** The values of DSCF at different positions of circular pipeline with *k*0*R* = 0.5 (*β*1 = 2.0).

**Table 2.** The values of DSCF at different positions of circular pipeline with *k*0*R* = 1.0 (*β*1 = 2.0).

**Table 3.** The values of DSCF at different positions of circular pipeline with *k*0*R* = 1.5 (*β*1 = 2.0).

In Figure 9, DSCF changes continuously with *β*1 at 30◦ of the circular pipeline with *k*0*R* = 1.0, 2.0 and 3.0. The *β*2 of Figure 9a–c are 0.5, 1.0 and 1.5, respectively. It can be seen that with the continuous increase of *β*1, the DSCF curves fluctuate significantly at the 30◦ of the circular pipeline. As the reference *k*0*R* increases, the oscillation frequency of DSCF curves increases. Another interesting point can be found in Figure 9b,c. Although the vibration frequencies of the DSCF curves in the subgraph are different, the maximum and minimum values in the same vibration cycle are approximately the same. It can be inferred that changing the density of the concrete with *β*1 can compensate for the influence of the reference wavenumber on the amplitude of DSCF. In addition, it can be observed that when the values of *β*2 increases, the oscillation frequency of DSCF curves will increase, but the oscillation amplitude of DSCF curves will decrease. It can be inferred that with the increase of concrete density and change of distribution, the amplitude of DSCF tends to be stable, and the fluctuation gradually decreases. The DSCF continuously changes curves with *β*1 is given in Figure 9d–f, *β*2 = −0.5, −1.0 and −1.5. SH waves incident horizontally from high density to low density of the concrete. It can be seen that compared with Figure 9a–c, the oscillation frequency of DSCF curves is relatively fast when the values of *β*2 are negative, and the minimum value of each vibration period is approaching 0. Similarly, the maximum and minimum values of the three DSCF curves under different reference wave numbers are approximately the same in the same vibration period.

In Figures 10 and 11, the DSCF curves continuously change with *β*2 at 30◦ of the circular pipeline *k*0*R* = 1.0, 2.0 and 3.0 are given. The values of *β*1 are 0.1, 0.3, 0.5. The values of *β*2 in Figure 10 are the positive of number of continuous changes with a range of 0.1~2.0, while *β*2 in Figure 11 are the negative number of a continuous range of −2.0~−0.1.

It can be observed in Figure 10 that the overall variation trend of DSCF curves increases with the increase of *β*2. When *β*1 = 0.1, it can be found that within the range of 0.1 < *β*2 < 1.5, the DSCF curves change with different values of *k*0*R* basically coincide without significant difference. When *β*2 > 0.5, the DSCF curves have obvious differences. When *β*1 = 0.3, the DSCF curves have obvious differences after *β*2 > 0.5. When *β*1 = 0.5, the DSCF curves have obvious differences after *β*2 > 0.2. Figure 11 shows the horizontal incident of SH waves from high density to low density of concrete. It can be found that as the absolute value of *β*2 increases, the DSCF curves at 30◦ of the circular pipeline have an obvious oscillation

phenomenon. When the values of *k*0*R* and *β*1 increase, the frequency of oscillation of DSCF curves will be higher, which is obviously different from the corresponding situation in Figure 10. When the values of *β*1 and *β*2 are both small, different *k*0*R* has little influence on the DSCF results at 30◦ of the circular pipeline. As the values of *β*1 increases the absolute value of *β*2, which causes differences between DSCF curves under different *k*0*R* conditions, will decrease, which is the same as the situation in Figure 10. This indicates that the density of the concrete with this density form will increase with the increases of *β*1, so that the DSCF curves with different wave numbers appear in advance.

**Figure 9.** The change of DSCF at the 30◦ position of the circular pipeline with the inhomogeneous parameter *β*1.

**Figure 10.** The change of DSCF at the 30◦ position of the circular pipeline with the inhomogeneous parameter *β*2 (0.1 ≤ *β*2 ≤ 2.0).

In Figure 12, the continuous change of DSCF curves with *k*0*R* at 30◦ of circular pipeline is given. In Figure 12a–c, *β*1 = 1.0, 1.5 and 2.0, *β*2 = 0.8, 1.0, 1.2. It can be observed that SH waves are incident horizontally from low density to high density of the concrete. It can be observed in Figure 12 that with the continuous change of *k*0*R*, the vibration amplitude of DSCF curves decreases. When the values of *β*1 and *β*2 increase, the oscillation frequency of DSCF curves will increase. In Figure 12d–f, the values of *β*2 = −0.8, −1.0 and −1.2, SH waves are incident horizontally from the high-density to the low-density direction of the concrete. In this case, compared with the case in Figure 12a–c, the maximum value of DSCF curves is less than, and the minimum points on the DSCF curves are all approaching 0. It can be found that similar to the continuous change of DSCF with *β*1 the maximum and minimum values of different DSCF curves in the same vibration period remain the same regardless of the values of *β*1 in each subgraph of Figure 12.

**Figure 11.** The change of DSCF at the 30◦ position of the circular pipeline with the inhomogeneous parameter *β*2 (−2.0 ≤ *β*2 ≤ −0.1).

In Figure 13, the continuous changes of DSCF curves with *k*0*R* at 30◦ of the circular pipeline are given. *β*2 = 1.2, 1.6, and 2.0, *β*1 = 0.1, 1.0 and 2.0. It can be observed that with the increase of *β*2, the peak value of DSCF also increases, the oscillation of DSCF curves will appear earlier and the oscillation frequency will be higher. At the same time, the larger *β*2 is, the closer the occurrence time of the maximum values of DSCF at the position of 30◦ is to the quasi-static condition (*k*0*R* = 0.1). With the continuous increase of *k*0*R*, the amplitude of DSCF curves decreases. In addition, larger *β*1 causes faster oscillation frequency of DSCF curves. Interestingly, even when the value of *β*1 changes, the DSCF curves with *β*2 = 1.2 always have an amplitude range of 4 to 6, the DSCF curves with *β*2 = 1.6 always have an amplitude range of 6 to 8, and the DSCF curves with *β*2 = 2.0 always have an amplitude range of 8 to 11.

**Figure 12.** The change of DSCF at the 30◦ position of the circular pipeline with the reference wave number *k*0*R* (*β*2 = ±0.8, ±1.0, ±1.2).

**Figure 13.** The change of DSCF at the 30◦ position of the circular pipeline with the reference wave number *k*0*R* (*β*1 = 0.1, 1.0, 2.0).

#### **7. Conclusions**

In this paper, the scattering of SH waves by circular pipeline in inhomogeneous concrete with polynomial-exponential coupling density distribution is studied based on the method of complex variable function. The analytical solution of dynamic stress concentration around circular pipeline is derived under this type of concrete with density variation. This paper discusses the effect of different dimensionless parameters on the distribution of DSCF around the circular pipeline. It provides a theoretical reference and a basis for analyzing the influence of defects on wave propagation in an inhomogeneous concrete with density variation. The specific conclusions are as follows:

(1) When the density inhomogeneous parameters *β*1, *β*2 and the reference wave number *k*0*R* change, all the maximum values of DSCF always appear on the back wave surface of the circular pipeline. In most cases, the maximum value of DSCF is concentrated within the range of 20~40◦ and 320~340◦ at the position of circular pipeline, and the number of extreme points of DSCF in this range also changes significantly.

The conformal transformation method used in this paper requires a high level of expression for the non-uniformity of concrete density, which is not yet achievable in reality. Although the inhomogeneous concrete model we presented does not exist in reality, we hope that in the future, the changing form and structural model of concrete density we proposed can be applied to concrete materials, and our research method can be applied to the elastic dynamic research of other inhomogeneous concrete forms.

**Author Contributions:** Conceptualization, Z.Y.; methodology, C.S.; software, G.J.; validation, Y.S. and X.L.; investigation, Y.Y.; data curation, C.S.; writing—original draft preparation, C.S.; writing—review and editing, C.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Scientific Research Fund of Institute of Engineering Mechanics, China Earthquake Administration, Grant No. 2021EEEVL0201, and the National Natural Science Foundation of China, Grant No U2239252 and the Research Team Project of Heilongjiang Natural Science Foundation, Grant No TD2020A001.

**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.

## **Appendix A. An Example of Choosing Values for Inhomogeneous Density Parameter and Reference Wave Number**

Reference [37] presents an inhomogeneous medium with a circular cavity, and its density expression is given as:

$$
\rho = \rho\_0 \left[ 4a^2 \left( \mathbf{x}^2 + y^2 \right) + 4a\beta \mathbf{x} + \beta^2 \right],
\tag{A1}
$$

where, *α* and *β* represent inhomogeneous parameters. When *β*/*α* < 5, the density of the medium is considered to vary continuously in two dimensions.

In the example, in order to analyze the variation of DSCF around a circular cavity under different reference wave numbers, reference wave numbers of 0.1, 0.5, 1.0, and 2.0 were selected. Therefore, in this article, we selected reference wave numbers of 0.5, 1.0, 1.5, and 2.0.

In the process of selecting values for inhomogeneous parameters in this example, it was controlled within the range of *β*/*α* < 5. The selected values for *α* and *β* are listed in Table A1.

Although the density variation form we provide is different from the density form in the example, we still give the range of variation for *β*1 values from 0.1 to 2.0 and *β*2 values from −2.0 to −0.1 and 0.1 to 2.0. However, it should be noted that the selection of values for inhomogeneous parameters can affect the accuracy of the results, and different values may be more suitable for different applications.

To illustrate the impact of inhomogeneous density and reference wave number values on the results, we analyzed the DSCF variation around a circular pipeline using different values for these parameters. We found that the selection of the density function exponent (*β*2 value) has a relatively large impact on the results compared to the values of *β*1 and reference wave number. This is because the selection of beta value can not only change the distribution form of density in the material, but also change the direction of density distribution.

In summary, our research results indicate the importance of carefully selecting the parameters used in inhomogeneous density concrete models, and emphasize the need for further research to better understand the impact of inhomogeneous density parameters and reference wave numbers on the DSCF obtained from these models.

**Table A1.** The values of *α* and *β*.

## **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* **Mixing Performance of a Passive Micro-Mixer with Mixing Units Stacked in Cross Flow Direction**

**Makhsuda Juraeva and Dong-Jin Kang \***

School of Mechanical Engineering, Yeungnam University, Gyoungsan 712-749, Korea<EMAIL_ADDRESS>**\*** Correspondence<EMAIL_ADDRESS>Tel.: +82-(0)53-810-2463

**Abstract:** A new passive micro-mixer with mixing units stacked in the cross flow direction was proposed, and its performance was evaluated numerically. The present micro-mixer consisted of eight mixing units. Each mixing unit had four baffles, and they were arranged alternatively in the cross flow and transverse direction. The mixing units were stacked in four different ways: one step, two step, four step, and eight step stacking. A numerical study was carried out for the Reynolds numbers from 0.5 to 50. The corresponding volume flow rate ranged from 6.33 µL/min to 633 µL/min. The mixing performance was analyzed in terms of the degree of mixing (*DOM*) and relative mixing energy cost (*MEC*). The numerical results showed a noticeable enhancement of the mixing performance compared with other micromixers. The mixing enhancement was achieved by two flow characteristics: baffle wall impingement by a stream of high concentration and swirl motion within the mixing unit. The baffle wall impingement by a stream of high concentration was observed throughout all Reynolds numbers. The swirl motion inside the mixing unit was observed in the cross flow direction, and became significant as the Reynolds number increased to larger than about five. The eight step stacking showed the best performance for Reynolds numbers larger than about two, while the two step stacking was better for Reynolds numbers less than about two.

**Keywords:** passive micro-mixer; mixing unit; cross flow direction; baffle impingement; swirl motion; mixing performance

## **1. Introduction**

Mixing at micro scales is an essential design task in the bio- and microfluidic systems such as micro-total analysis system (µ-TAS), lab on a chip, and micro-reactors. As the associated dimensions of a microfluidic system are small, the molecular diffusion is a major mechanism of mixing; the corresponding Reynolds number is also very small and the flow is laminar. Meanwhile, many biochemical processes require rapid and complete mixing, and mixing enhancement is essential in designing these systems [1].

A variety of micro-mixers have been devised to enhance the mixing in a micro-fluidic system [2]. They are named either active or passive micro-mixers, depending on the usage of an external energy source. As active micro-mixers utilize an external energy source to achieve the mixing enhancement, the structure of a micro-mixer becomes more complicated than that of passive micro-mixers. In addition, it is more costly compared with a passive micro-mixer. Examples of external energy sources used for active micro-mixers are acoustic [3], magnetic [4], electric [5], thermal [6], electro [7,8], and pressure fluctuating [9]. On the contrary, passive micro-mixers use micro-mixer geometry to agitate and generate secondary flow. They are expected to promote chaotic advection of fluids leading to the mixing enhancement. Accordingly, the mixing enhancement by a passive micro-mixer could be limited compared with that by a similar, active micro-mixer. However, passive micro-mixers are simpler to fabricate and easier to integrate into microfluidic systems; therefore, they are widely used in microfluidic systems.

Many different types of passive micro-mixers have been designed and investigated to enhance their mixing performance [10]. Some examples include recessed grooves in the

**Citation:** Juraeva, M.; Kang, D.-J. Mixing Performance of a Passive Micro-Mixer with Mixing Units Stacked in Cross Flow Direction. *Micromachines* **2021**, *12*, 1530. https://doi.org/10.3390/mi12121530

Academic Editor: Kwang-Yong Kim

Received: 25 October 2021 Accepted: 6 December 2021 Published: 9 December 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/).

channel wall [11], herringbone-type walls [12], baffles [13,14], channel wall twisting [15], junction contraction [16], split and recombine (SAR) [17,18], serpentine microchannel [19,20], converging and diverging microchannel [21], spiral channel [22,23], and periodic geometric features [24].

Among various geometric components, baffles are one of the most widely adopted components in designing a passive micro-mixer. Many researchers have shown that micromixers based on baffles perform effectively over a wide range of Reynolds numbers [13,25]. For example, Tsai et al. [25] used radial baffles in a curved microchannel and showed that radial baffles induced vortices in multiple directions. Kang [13] placed rectangular baffles in a cyclic order along the channel wall and enhanced the mixing performance of a Tshaped micromixer. Santana et al. [14] used triangular baffles to obtain a high mixing index for various micro-channel dimensions. Sotowa et al. [26] showed that the indentations and baffles attached to the micromixer wall enhance the mixing by means of secondary flow in deep micro-channel reactors. Chen et al. [27] distributed baffles on both side walls of a micromixer based on the Koch fractal principle. They showed that the baffles distributed on the both side of the microchannel wall performs better than the mixing performance of the baffles distributed on the one side. Chung et al. [28] implemented planar baffles and showed that the mixing performance was enhanced at both of diffusion dominant and convection dominant regime of mixing. Raza et al. [29] improved the mixing performance of a SAR micromixer by embedding baffles just after each SAR unit. The enhancement was shown as achievable over a Reynolds number range from 0.1 to 80.

Recently, several complicated three-dimensional (3D) micro-mixers showed an improved mixing performance compared with those of two-dimensional (2D) micro-mixers of similar size [30–33]. However, 2D planar micromixers have an advantage of simplicity in fabrication compared with that of complex 3D micromixers. Various design techniques were attempted to promote 3D effects: baffles, channel wall twisting, spiral channel, and split and recombine. For example, Kang [13] showed that a cyclic configuration of baffles generates a rotational flow in the cross section of a micromixer. The same idea was adopted in the present micromixer to promote 3D effects, and four baffles were attached to the four side walls of the micromixer in a cyclic order. The overall layout of a passive micromixer is also an important design concept [33]. For example, Tripathi [23] studied three different layout of a passive micromixer based on spiral microchannel and showed that the mixing performance is quite dependent on the micromixer layout. The eight mixing units were stacked in the cross flow direction; this kind of stacking has not been studied yet. The stacking was designed as four different layouts: one step, two step, four step, and eight step stacking. The combined effects of a cyclic baffle arrangement and mixing unit stacking in the cross flow direction was the main design concept of the present micromixer.

Most of micromixers used in biological and chemical applications usually operate in the range of milliseconds of mixing time, and the corresponding Reynolds number is less than about 100 [34–36]; the corresponding volume flow rate is an order of mL/h. In this range, the mixing mechanism can be divided into three regimes: molecular dominance, transition, and advection dominance [10,18,33]. The effects of micromixer layout design as well as a passive device such as baffles are known to be significant in all of the three regimes. Accordingly, a numerical study was carried out for the Reynolds numbers from 0.5 to 50, covering all of the three regimes. The corresponding volume flow rate ranged from 6.33 µL/min to 633 µL/min.

A numerical simulation has several advantages to study the fluid dynamic and mixing features involved in a micro-mixer, including easy visualization of the mixing process and the associated flow characteristics such as streamlines and vortex formation. Accordingly, a numerical approach is widely accepted in studying mixing enhancement mechanism of a micro-mixer. For example, Rhoades et al. [37] used the commercial software COMSOL Multiphysics 5.1 (COMSOL, Inc., Burlington, MA, USA) to study the mixing performance of grooved serpentine microchannels. Volpe et al. [38] studied the flow dynamics of a continuous size-based sorter microfluidic device by using the lattice Boltzmann method (LBM).

Kang [13] used the commercial software ANSYS® Fluent 19.2 (ANSYS, Inc., Canonsburg, PS, USA) [39] to simulate the mixing performance of a passive micromixer with baffles and quantitatively evaluated the mixing performance. In this paper, the commercial software ANSYS® Fluent 19.2 was used to simulate the mixing performance. The mixing performance was estimated in terms of the degree of mixing (*DOM*) and corresponding mixing energy cost (*MEC*). The numerical simulation

micromixer with baffles and quantitatively evaluated the mixing performance.

and the associated flow characteristics such as streamlines and vortex formation. Accordingly, a numerical approach is widely accepted in studying mixing enhancement mechanism of a micro-mixer. For example, Rhoades et al. [37] used the commercial software COMSOL Multiphysics 5.1 (COMSOL, Inc., Burlington, MA, USA) to study the mixing performance of grooved serpentine microchannels. Volpe et al. [38] studied the flow dynamics of a continuous size-based sorter microfluidic device by using the lattice Boltzmann method (LBM). Kang [13] used the commercial software ANSYS® Fluent 19.2 (AN-SYS, Inc., Canonsburg, PS, USA) [39] to simulate the mixing performance of a passive

*Micromachines* **2021**, *12*, x FOR PEER REVIEW 3 of 18

In this paper, the commercial software ANSYS® Fluent 19.2 was used to simulate the mixing performance. The mixing performance was estimated in terms of the degree of mixing (*DOM*) and corresponding mixing energy cost (*MEC*). The numerical simulation was carried out for Reynolds numbers ranging from 0.5 to 50, corresponding to volume flow rates ranging from 6.33 µL/min to 633 µL/min. was carried out for Reynolds numbers ranging from 0.5 to 50, corresponding to volume flow rates ranging from 6.33 μL/min to 633 μL/min. **2. Passive Micro-Mixer with Mixing Units Stacked in the Cross Flow Direction**  Figure 1 shows a mixing unit with four hexahedron baffles inside. The four baffles

#### **2. Passive Micro-Mixer with Mixing Units Stacked in the Cross Flow Direction** were arranged in a cyclic order. The first baffle was attached to the lower wall, and the

Figure 1 shows a mixing unit with four hexahedron baffles inside. The four baffles were arranged in a cyclic order. The first baffle was attached to the lower wall, and the third baffle was attached to the upper wall. The second and fourth baffles were attached to the front and back walls, respectively. Each baffle was 30 µm thick, and blocked the flow passage by half in each direction. Two consecutive baffles were separated by 60 µm, and each mixing unit was 600 µm long. Figure 2 depicts the five different layouts of a micromixer simulated in this paper. Each micromixer had eight mixing units, and each mixing unit consisted of four baffles. Figure 2a shows the baseline layout of eight mixing units, and Figure 2b–e shows mixing units stacked in the cross stream direction (y direction); the mixing units were overlapped 150 µm in the main stream direction (x direction). The inlets and outlet had a rectangular cross-section of 300 µm × 120 µm. The two inlet branches were 1000 µm long, and the outlet branch was 600 µm long. third baffle was attached to the upper wall. The second and fourth baffles were attached to the front and back walls, respectively. Each baffle was 30 μm thick, and blocked the flow passage by half in each direction. Two consecutive baffles were separated by 60 μm, and each mixing unit was 600 μm long. Figure 2 depicts the five different layouts of a micromixer simulated in this paper. Each micromixer had eight mixing units, and each mixing unit consisted of four baffles. Figure 2a shows the baseline layout of eight mixing units, and Figure 2b–e shows mixing units stacked in the cross stream direction (y direction); the mixing units were overlapped 150 μm in the main stream direction (x direction). The inlets and outlet had a rectangular cross-section of 300 μm × 120 μm. The two inlet branches were 1000 μm long, and the outlet branch was 600 μm long. For the sake of simplicity, we assumed that the same aqueous solution flows into the

For the sake of simplicity, we assumed that the same aqueous solution flows into the two inlets: Inlet 1 and Inlet 2. The properties of the aqueous solution are similar to those used in BioMEMS systems. The kinematic viscosity and mass diffusion coefficient of the aqueous solution are *v* = 10−6 m2/s and *D* = 10−9 m2/s, respectively [18,33]. The corresponding Schmidt (Sc) number is 103 ; it is defined as the ratio of the kinetic viscosity and the mass diffusion coefficient. The Reynolds number is defined as *Re* = *Umeandh v* , where *Umean*, *dh* and *v* indicate the mean velocity at the outlet, the hydraulic diameter of outlet, and the kinematic viscosity of the fluid, respectively. two inlets: Inlet 1 and Inlet 2. The properties of the aqueous solution are similar to those used in BioMEMS systems. The kinematic viscosity and mass diffusion coefficient of the aqueous solution are = 10−6 m2/s and *D* = 10−9 m2/s, respectively [18,33]. The corresponding Schmidt (Sc) number is 103; it is defined as the ratio of the kinetic viscosity and the mass diffusion coefficient. The Reynolds number is defined as = ೌௗ ௩ , where , and indicate the mean velocity at the outlet, the hydraulic diameter of outlet, and the kinematic viscosity of the fluid, respectively.

**Figure 1. Figure 1.**  Schematic diagram of the four-baffle arrangement within a mixing unit. Schematic diagram of the four-baffle arrangement within a mixing unit.

**Figure 2.** Present passive micro-mixers with mixing units stacked in the cross flow direction. (**a**) Baseline layout. (**b**) One step stacking. (**c**) Two step stacking. (**d**) Four step stacking. (**e**) Eight step stacking. **Figure 2.** Present passive micro-mixers with mixing units stacked in the cross flow direction. (**a**) Baseline layout. (**b**) One step stacking. (**c**) Two step stacking. (**d**) Four step stacking. (**e**) Eight step stacking.

#### **3. Governing Equations and Computation Procedure 3. Governing Equations and Computation Procedure**

Transport phenomena in micromixers can be described theoretically at two different levels: the molecular level and the continuum level. The two different levels are related to the typical length scale involved. The continuum model can describe most transport phenomena in micromixers with a length scale ranging from micrometers to centimeters [40]. The present micromixers were in this range of length scale. Transport phenomena in micromixers can be described theoretically at two differentlevels: the molecular level and the continuum level. The two different levels are related tothe typical length scale involved. The continuum model can describe most transport phenomena in micromixers with a length scale ranging from micrometers to centimeters [40]. The present micromixers were in this range of length scale.

In this study, the fluid entering the two inlets were water and dyed water. They were assumed to have the same physical properties of water at a temperature of 20 °C; the fluid A at Inlet 1 was the dyed water and the fluid B at Inlet 2 was water. The flow inside the micromixer was assumed to be steady, incompressible, laminar, and Newtonian. There-In this study, the fluid entering the two inlets were water and dyed water. They were assumed to have the same physical properties of water at a temperature of 20 ◦C; the fluid A at Inlet 1 was the dyed water and the fluid B at Inlet 2 was water. The flow inside the micromixer was assumed to be steady, incompressible, laminar, and Newtonian. Therefore,

fore, the fluid flow in the micromixer was simulated solving the incompressible Navier–

Stokes and the continuity equations:

the fluid flow in the micromixer was simulated solving the incompressible Navier–Stokes and the continuity equations:

$$\left(\overrightarrow{\boldsymbol{\mu}} \cdot \nabla\right)\overrightarrow{\boldsymbol{\mu}} = -\frac{1}{\rho} \nabla p + \mathbf{v} \nabla^2 \overrightarrow{\boldsymbol{\mu}}\tag{1}$$

$$\nabla \cdot \stackrel{\rightarrow}{u} = \mathbf{0} \tag{2}$$

where *ρ*, → *u*, *p*, and ν are the density, the velocity vector, the pressure, and the kinematic viscosity of the fluid, respectively. The density and viscosity were specified as 998 kg/m3 and 0.001 kg/(ms), respectively.

The flow field thus obtained was used to simulate the mixing process throughout the micromixer. In this study, the mixing was assumed to be governed by two fluid dynamic mechanisms: molecular diffusion and advection. A scalar advection-diffusion transport equation was used to simulate the mixing process [33,37]:

$$\left(\stackrel{\rightarrow}{u}\nabla\right)\phi = D\nabla^2\phi\tag{3}$$

where *D* and *φ* are the diffusion coefficient and the dye concentration, respectively. The dyed water concentration is *φ* = 1 at Inlet 1, and *φ* = 0 at Inlet 2. In general, a concentration gradient is a density gradient which can induce a convective flow. The flow velocity due to concentration gradients can vary from 0.1 to 10 µms−1 depending on properties such as the channel dimensions, fluid viscosity, fluid type, gradient magnitude [41]; in this paper, the flow velocity in the micromixer was in the range of 0.001 to 0.14 ms−1 . Therefore, the convective flow induced by any concentration gradient was neglected in this paper.

We used the commercial software ANSYS® Fluent to solve the governing Equations (1)–(3), and it is based on the finite volume method. In general, a certain amount of numerical diffusion is introduced in discretizing the convective terms, but it can be limited by using a high-order discretization scheme. The non-linear convective terms in Equations (1) and (3) were approximated using the QUICK scheme (quadratic upstream interpolation for convective kinematics), and its theoretical accuracy is third order. The velocity was assumed uniform at the two inlets, while it was assumed fully developed at the outlet. Along all the other walls, the no-slip boundary was applied.

In order to obtain a fully converged solution, every computation was continued until oscillation of the residual of all equations was negligibly small. Some researchers used an adaptive mesh technique to reduce the oscillation [42]. In this study, the residuals of continuity, momentum, and concentration equations oscillated in the range smaller than 10−11, 10−14, and 10−8 , respectively; they were sufficiently small to obtain reliable numerical solutions.

Accurate numerical simulation of the mixing in micromixers is still a challenging problem, especially for high Peclet numbers. Some studies do not deal with computational issues related to the mesh dependence of numerical solution. In this paper, a detailed study was conducted in the section of validation of numerical solution. According to Okuducu et al. [43], the accuracy of numerical solutions is also dependent on the type of cells; structured hexahedral cells show the most reliable numerical solution, in comparison with tetrahedral and prism cells. In this paper, all cells for every micromixer layout were structured and hexahedral.

The mixing performance of present micro-mixer was evaluated using the degree of mixing (*DOM*) and mixing energy cost (*MEC*). We calculated the *DOM* in the following form:

$$DOM = 1 - \frac{1}{\tilde{\mathfrak{F}}} \sqrt{\sum\_{i=1}^{n} \frac{\left(\phi\_i - \tilde{\mathfrak{F}}\right)^2}{n}} \tag{4}$$

where *ϕι* and *n* are the mass fraction of fluid A in the *i*th cell and total number of cells, respectively, and *ξ* is 0.5 when the two fluids are completely mixed. The *MEC* is calculated in the following form, and measures the effectiveness of a micro-mixer [44,45]:

∆

*Micromachines* **2021**, *12*, x FOR PEER REVIEW 6 of 18

$$\text{MEC} = \frac{\Delta p / \rho u\_{\text{mean}}^2}{DOMX100} \tag{5}$$

where *umean* is the average velocity at the outlet, and ∆*p* is the pressure difference between the inlet and the outlet. where is the average velocity at the outlet, and ∆ is the pressure difference between the inlet and the outlet.

#### **4. Validation of Numerical Study 4. Validation of Numerical Study**

For high Sc number simulations, numerical diffusion is known to deteriorate the accuracy of simulated results, in general [46–49]. To minimize the numerical diffusion problem, several approaches can be chosen. These include a particle-based simulation such as the Monte Carlo method [46] or a grid-based method with a small cell Peclet number. Here, the cell Peclet number is *Pe* = *Ucell lcell D* where *Ucell* and *lcell* are the local flow velocity and cell size, respectively. However, these approaches are computationally too expensive to adopt in a study like this paper. Most numerical studies prefer a practical approach to obtain numerical solutions with a reasonable degree of accuracy. To achieve that, a detailed study of grid independence including the grid convergence index (GCI) test is usually adopted [18,29,31]. In this paper, a similar procedure was followed. For high Sc number simulations, numerical diffusion is known to deteriorate the accuracy of simulated results, in general [46–49]. To minimize the numerical diffusion problem, several approaches can be chosen. These include a particle-based simulation such as the Monte Carlo method [46] or a grid-based method with a small cell Peclet number. Here, the cell Peclet number is = where and are the local flow velocity and cell size, respectively. However, these approaches are computationally too expensive to adopt in a study like this paper. Most numerical studies prefer a practical approach to obtain numerical solutions with a reasonable degree of accuracy. To achieve that, a detailed study of grid independence including the grid convergence index (GCI) test is usually adopted [18,29,31]. In this paper, a similar procedure was followed.

To quantitatively validate the present numerical approach, we simulated a passive micro-mixer experimented by Tsai et al. [25]. Figure 3 shows a schematic diagram of the micro-mixer. The two inlets had a rectangular cross section of width 45 µm by depth 130 µm. The micromixer had four baffles of width 45 µm by height 97.5 µm. The fluid was assumed to have the properties of density 997 kg/m3 , viscosity 0.00097 kg/(ms), and diffusion coefficient 3.6 × 10−10 m2/s, as reported by Tsai et al. [25]. The numerical simulation was carried out for three different Reynolds numbers of *Re =* 1, 9, and 81 and the results were compared with the corresponding experimental data. To quantitatively validate the present numerical approach, we simulated a passive micro-mixer experimented by Tsai et al. [25]. Figure 3 shows a schematic diagram of the micro-mixer. The two inlets had a rectangular cross section of width 45 μm by depth 130 μm. The micromixer had four baffles of width 45 μm by height 97.5 μm. The fluid was assumed to have the properties of density 997 kg/m3, viscosity 0.00097 kg/(ms), and diffusion coefficient 3.6 × 10−10 m2/s, as reported by Tsai et al. [25]. The numerical simulation was carried out for three different Reynolds numbers of *Re =* 1, 9, and 81 and the results were compared with the corresponding experimental data.

**Figure 3.** Diagram of the micromixer experimented by Tsai et al. [25]. **Figure 3.** Diagram of the micromixer experimented by Tsai et al. [25].

Hexahedral cells were used to mesh the computational domain sketched in Figure 3; they are all structured. Before detailed simulations, a preliminary study was carried out to check the grid independence of numerical solutions for the Reynolds number of 9. Figure 4a shows the dependence of the numerical solution for *Re* = 1, and about 1 million of the computational cells were enough to obtain a numerical solution with a reasonable accuracy. Here, *DOMT* stands for the degree of mixing defined by Tsai et al. [25] in the following way: Hexahedral cells were used to mesh the computational domain sketched in Figure 3; they are all structured. Before detailed simulations, a preliminary study was carried out to check the grid independence of numerical solutions for the Reynolds number of 9. Figure 4a shows the dependence of the numerical solution for *Re* = 1, and about 1 million of the computational cells were enough to obtain a numerical solution with a reasonable accuracy. Here, *DOMT* stands for the degree of mixing defined by Tsai et al. [25] in the following way:

$$DOM\_T = 1 - \frac{\sigma\_D}{\sigma\_{D,\rho}} \tag{6}$$

and

and by Tsai et al. [25] for Reynolds numbers from 1 to 81. The numerical solution and experi-

and
$$\sigma = \sqrt{\frac{1}{n} \sum\_{i=1}^{n} (\phi\_i - \phi\_{ave})^2} \tag{7}$$

(7)

*Micromachines* **2021**, *12*, x FOR PEER REVIEW 7 of 18

*Micromachines* **2021**, *12*, x FOR PEER REVIEW 7 of 18

= ඩ1

where *σD* is the standard deviation of *φ* on a cross section normal to the flow, *σD*,*o* is the standard deviation at the inlet, and *φave* is the average value of *φ* at a sampled cross section. the Peclet number decreased. The discrepancy is attributed to several factors such as the numerical diffusion, experimental uncertainty, etc. where is the standard deviation of on a cross section normal to the flow, , is the standard deviation at the inlet, and ௩ is the average value of at a sampled cross section. Figure 4b compares the simulation results with the corresponding experimental data

ሺ − ௩ሻଶ

ୀଵ

where is the standard deviation of on a cross section normal to the flow, , is the standard deviation at the inlet, and ௩ is the average value of at a sampled cross section. Figure 4b compares the simulation results with the corresponding experimental data

**Figure 4.** Validation of numerical solution. (**a**) Grid dependence of the numerical solution. (**b**) *DOM* vs. *Re*. **Figure 4.** Validation of numerical solution. (**a**) Grid dependence of the numerical solution. (**b**) *DOM* vs. *Re*.

Prior to carrying out the final simulations, an additional set of preliminary simulations was carried out to determine an appropriate mesh size for the simulation of the present micro-mixer. We used the present micro-mixer of Figure 2d for this study, and simulated for the Reynolds number of *Re* = 3. The size of every hexahedral cell was varied from 2.5 μm to 5 μm. Figure 5 shows the dependence of the calculated *DOM* on the edge size. The deviation of 3.5 μm solution from that of 4 μm was 0.4%. Therefore, 3.5 μm was small enough to obtain grid independent solutions. Figure 4b compares the simulation results with the corresponding experimental data by Tsai et al. [25] for Reynolds numbers from 1 to 81. The numerical solution and experimental data showed similar behavior as the Reynolds number increased, and the discrepancy was acceptable. The discrepancy between the experimental data and numerical solution was less than 4%, and became smaller as the Reynolds number decreased, and as the Peclet number decreased. The discrepancy is attributed to several factors such as the numerical diffusion, experimental uncertainty, etc. (**a**) (**b**) **Figure 4.** Validation of numerical solution. (**a**) Grid dependence of the numerical solution. (**b**) *DOM* vs. *Re*.

Prior to carrying out the final simulations, an additional set of preliminary simulations was carried out to determine an appropriate mesh size for the simulation of the present micro-mixer. We used the present micro-mixer of Figure 2d for this study, and simulated for the Reynolds number of *Re* = 3. The size of every hexahedral cell was varied from 2.5 µm to 5 µm. Figure 5 shows the dependence of the calculated *DOM* on the edge size. The deviation of 3.5 µm solution from that of 4 µm was 0.4%. Therefore, 3.5 µm was small enough to obtain grid independent solutions. Prior to carrying out the final simulations, an additional set of preliminary simulations was carried out to determine an appropriate mesh size for the simulation of the present micro-mixer. We used the present micro-mixer of Figure 2d for this study, and simulated for the Reynolds number of *Re* = 3. The size of every hexahedral cell was varied from 2.5 μm to 5 μm. Figure 5 shows the dependence of the calculated *DOM* on the edge size. The deviation of 3.5 μm solution from that of 4 μm was 0.4%. Therefore, 3.5 μm was small enough to obtain grid independent solutions.

**Figure 5. Figure 5.** Dependence of the numerical solution on the edge size. Dependence of the numerical solution on the edge size.

Using the preliminary simulation results, the uncertainty of grid convergence was evaluated using the grid convergence index (*GCI*) [50,51]. According to the Richardson extrapolation methodology, the *GCI* is defined as follows:

$$GCI = F\_s \frac{|\varepsilon|}{r^p - 1} \tag{8}$$

$$
\varepsilon = \frac{f\_{\text{coarse}\varepsilon} - f\_{\text{fine}}}{f\_{\text{fine}}} \tag{9}
$$

where *Fs* , *r*, and *p* are the safety factor of the method, grid refinement ratio, and the order of accuracy of the numerical method, respectively. *fcoarse* and *ffine* are the numerical results obtained with a coarse grid and fine grid, respectively. We specified *Fs* as 1.25 according to the suggestion by Roache [50]. For the GCI evaluation, the edge size was varied from 2.5 µm to 5 µm, 2.5 µm, 3.5 µm, 4 µm, and 5 µm. Table 1 summarizes the result of the *GCI* test. The *GCI* based on the *DOM* at the outlet was 1% and 0.57% for 4 µm and 3.5 µm, respectively. Therefore, the edge size of 3.5 µm was small enough to obtain numerical solutions with reasonable accuracy.

**Table 1.** Results of the GCI test.

#### **5. Results**

A new passive micro-mixer with mixing units stacked in the cross flow direction was proposed, its mixing performance was simulated for Reynolds numbers ranging from 0.1 to 50. The present micro-mixer consisted of eight mixing units, and each mixing unit contained four baffles; the baffles were arranged in three different ways. The mixing units were stacked in four different ways to study their effect on the mixing performance.

The velocity at the two inlets was uniform in the range from 0.0293 mm/s to 146 mm/s, and the corresponding volume flow rates were from 6.33 to 633 µL/min. The mixing performance was evaluated in terms of the *DOM* at the outlet and the corresponding *MEC*.

Figure 6 shows the effects of mixing unit stacking in the cross flow direction on the mixing performance in terms of the *DOM* and the corresponding *MEC*. All of the cross flow stacking layouts resulted in a noticeable variation from the baseline layout, and the variation was prominent in the Reynolds number range from 1 to 30. For the Reynolds numbers less than about one, the two step stacking showed the best performance in terms of *DOM*. For example, the *DOM* of the two step stacking at *Re* = 1 showed a 9% increase from that of the baseline layout. On the contrary, the eight step stacking showed the best *DOM* for the Reynolds numbers larger than about two. For example, the *DOM* of the eight step stacking at *Re* = 20 was enhanced 14% from that of the baseline layout. However, the eight step stacking showed the least effective performance for the Reynolds numbers less than about one, as can be seen in the *MEC* distribution; it required the largest pressure load.

**Figure 6.** Effects of mixing unit stacking on the mixing performance. (**a**) *DOM* vs. Re. (**b**) *MEC* vs. Re. **Figure 6.** Effects of mixing unit stacking on the mixing performance. (**a**) *DOM* vs. *Re*. (**b**) *MEC* vs. *Re*.

We examined how the mixing unit stacking affected the mixing performance of a passive micro-mixer. Figure 7 shows the distribution of mass concentration of the fluid A at the plane of channel half width *z* = 60 μm for *Re* = 1. The concentration contours showed that a stream of high concentration of fluid A (reddish yellow line in the figure) formed along the center of the flow passage. Its appearance from the second mixing unit suggests that the mixing in the first mixing unit was actively achieved by the baffles. One interesting thing to note is that the stream of high concentration impinged on the baffle in the mixing units, especially in the two step stacking layout. Considering that the convective mixing effects for *Re* = 1 were quite limited, the flow characteristics of baffle impingement is an interesting phenomenon in the diffusion dominant flow regime. The baffle impingement seemed most intensive for the case of two step stacking. To quantitatively check this explanation, the *DOM* increment in each mixing unit was compared for the baseline and two step stacking layouts. Figure 8 shows the *DOM* increment in each mixing unit for *Re*  = 1. Here, the *DOM* increment means the difference of *DOM* between at the exit and entry of each mixing unit. For example, the flow entry and exit in the second mixing unit of the two step layout are indicated as the arrows pointing upward and downward (see Figure 2c), respectively. The entry and exit *DOMs* were calculated at the corresponding cross section normal to the flow. The increment of the two step stacking in the third and fourth mixing units was noticeably larger than that of the baseline layout. For example, the *DOM* increment in the third mixing unit of two step stacking was about 64% larger than the baseline layout. We examined how the mixing unit stacking affected the mixing performance of a passive micro-mixer. Figure 7 shows the distribution of mass concentration of the fluid Aat the plane of channel half width *z* = 60 µm for *Re* = 1. The concentration contours showed that a stream of high concentration of fluid A (reddish yellow line in the figure) formed along the center of the flow passage. Its appearance from the second mixing unit suggests that the mixing in the first mixing unit was actively achieved by the baffles. One interesting thing to note is that the stream of high concentration impinged on the baffle in the mixing units, especially in the two step stacking layout. Considering that the convective mixing effects for *Re* = 1 were quite limited, the flow characteristics of baffle impingement is an interesting phenomenon in the diffusion dominant flow regime. The baffle impingement seemed most intensive for the case of two step stacking. To quantitatively check this explanation, the *DOM* increment in each mixing unit was compared for the baseline and two step stacking layouts. Figure 8 shows the *DOM* increment in each mixing unit for *Re* = 1. Here, the *DOM* increment means the difference of *DOM* between at the exit and entry of each mixing unit. For example, the flow entry and exit in the second mixing unit of the two step layout are indicated as the arrows pointing upward and downward (see Figure 2c), respectively. The entry and exit *DOMs* were calculated at the corresponding cross section normal to the flow. The increment of the two step stacking in the third and fourth mixing units was noticeably larger than that of the baseline layout. For example, the *DOM* increment in the third mixing unit of two step stacking was about 64% larger than the baseline layout.

Figure 9 shows the distribution of mass concentration of the fluid A at the plane of channel half width *z* = 60 μm for *Re* = 20. The concentration contours showed a very chaotic distribution in the first and second mixing units. However, the stream of fluid A (red stream in the figure) and fluid B (blue stream in the figure) started to segregate from each other in the third mixing unit for the baseline layout. On the other hand, the eight step stacking layout showed a quite different flow pattern. The stream of fluid B impinged on the baffle in the third and fourth mixing units while the stream of fluid A impinged on the baffle in the second, fifth, and sixth mixing units. This suggests that the stream of high concentration changes their position as it flows down along the mixing unit. Figure 10 plots the distribution of mass concentration of the fluid A at the cross sections before and after each mixing unit for the eight step stacking layout, with contours in the *zx* plane (in the cross flow direction). The streams of fluid A and B showed a swirl motion as they passed through the mixing units. Since the swirl motion was in the y-direction (*zx* plane), the stacking of the mixing units in the cross flow direction contributed to the mixing enhancement discussed in Figure 6. This suggests that the mixing enhancement in the convection dominant flow regime can be easily obtained by a simple stacking method of the Figure 9 shows the distribution of mass concentration of the fluid A at the plane of channel half width *z* = 60 µm for *Re* = 20. The concentration contours showed a very chaotic distribution in the first and second mixing units. However, the stream of fluid A (red stream in the figure) and fluid B (blue stream in the figure) started to segregate from each other in the third mixing unit for the baseline layout. On the other hand, the eight step stacking layout showed a quite different flow pattern. The stream of fluid B impinged on the baffle in the third and fourth mixing units while the stream of fluid A impinged on the baffle in the second, fifth, and sixth mixing units. This suggests that the stream of high concentration changes their position as it flows down along the mixing unit. Figure 10 plots the distribution of mass concentration of the fluid A at the cross sections before and after each mixing unit for the eight step stacking layout, with contours in the *zx* plane (in the cross flow direction). The streams of fluid A and B showed a swirl motion as they passed through the mixing units. Since the swirl motion was in the y-direction (*zx* plane), the stacking of the mixing units in the cross flow direction contributed to the mixing enhancement discussed in Figure 6. This suggests that the mixing enhancement in the convection dominant flow regime can be easily obtained by a simple stacking method of the mixing units in the cross flow direction.

mixing units in the cross flow direction.

**Figure 7. Figure 7.**  Concentration distribution for Concentration distribution for *Re* = 1 at the plane of half width *Re* = 1 at the plane of half width *z* = 60*z*µ = 60 m. μm.

*Micromachines* **2021**, *12*, x FOR PEER REVIEW 11 of 18

**Figure 8.** Comparison of *DOM* increment in each mixing unit for *Re* = 1. **Figure 8.** Comparison of *DOM* increment in each mixing unit for *Re* = 1.

Figure 11 shows the velocity vector field at the junction of third and fourth mixing

**Figure 9. Figure 9.**  Concentration distribution for Concentration distribution for *Re Re* = 20 at the plane of half width = 20 at the plane of half width *z* = 60 μm. *z* = 60 µm.

**Figure 10.** Concentration distribution at the cross section before and after the mixing unit for *Re* = 20. **Figure 10.** Concentration distribution at the cross section before and after the mixing unit for *Re* = 20.

**Figure 11.** Velocity vector field at the junction of the third and fourth mixing units. Figure 13 compares the simulated mixing performance in terms of the *DOM* and re-Figure 11 shows the velocity vector field at the junction of third and fourth mixing units (at section (d) in Figure 10) for *Re* = 20, 5, and 2. The velocity vector field in Figure 11a shows two counter rotating vortices. A larger and stronger vortex flow was observed on the lower right corner, and the other vortex was on the upper left corner. Comparing this with Figure 10d, the vortex flow pattern caused a stream of high concentration of fluid A (a red stream on the lower right corner) swirling in the clockwise direction along the mixing unit walls. However, this swirl motion became weak as the Reynolds number decreased, as can be seen in Figure 11b,c. A weak vortex seemed to form near the lower right corner for *Re* = 5, while no explicit pattern of vortex was observed for *Re* = 2. This suggests that the mixing unit stacking in the cross flow direction generated a swirl motion inside the mixing unit, even if the swirl motion was insignificant for Reynolds numbers less than about five. The swirl motion inside the mixing unit is the second flow characteristic attributed to the mixing enhancement described in the Figure 8. **Figure 10.** Concentration distribution at the cross section before and after the mixing unit for *Re* = 20.

bination of the cross channel split and recombine (CC-SAR) and the mixing cell with baffles **Figure 11.** Velocity vector field at the junction of the third and fourth mixing units. **Figure 11.** Velocity vector field at the junction of the third and fourth mixing units.

(MC-B). They showed that the mixing performance of the optimized micromixer was noticeably improved compared with a micromixer based on the CC-SAR alone. Hossain et al. [53] studied a micromixer with unbalanced three split rhombic sub-channels. According to their Figure 13 compares the simulated mixing performance in terms of the *DOM* and required pressure load with those from other passive micromixers. In general, the mixing performance of micromixers increased with the diffusion coefficient. The value of diffusion coefficient depends on the fluid used in actual applications. To neglect the effects of the diffusion coefficient on the mixing performance, several passive micromixers based on the same size of Figure 12 compares the streamlines starting from the two inlets and concentration contours at a plane within each mixing unit for the baseline and eight step stacking layouts for *Re* = 20. The streamline plots were obtained at the cross section of half width *z* = 60 µm. On the right upper side, an enlarged view of the streamlines passing through from the third to fifth mixing units are shown. The concentration contours were obtained at the

the diffusion coefficient were selected for comparison; the diffusion coefficient was 1.2 × 10−9 m2/s. Tripathi et al. [23] numerically investigated the mixing performance of spiral micromix-

tices was found to enhance the mixing performance in the spiral micromixer. Raza et al. [29] proposed a SAR mixing unit combined with baffles in a curved channel and simulated its mixing performance. They showed that the mixing performance was better than their earlier version of SAR micro-mixer, and attributed the mixing enhancement to the combined effects of collision, Dean vortex, and secondary flows. Makhsuda et al. [52] proposed an optimal combination of the cross channel split and recombine (CC-SAR) and the mixing cell with baffles (MC-B). They showed that the mixing performance of the optimized micromixer was noticeably improved compared with a micromixer based on the CC-SAR alone. Hossain et al. [53] studied a micromixer with unbalanced three split rhombic sub-channels. According to their plane of the first baffle within each mixing unit along the micromixer. The streamlines from the inlet 1 (red) and streamlines from the inlet 2 (blue) changed their relative position more frequently for the eight step stacking layout as they passed down the mixing units; this flow pattern was more clearly seen in the enlarged view shown on the right upper side. The rapid change of their relative position led to faster mixing at the cross section, as can be seen in the centration contours within the mixing units. This flow pattern is also explained by the velocity vector field shown on the right lower side. The eight step stacking layout showed a vortex flow on the right lower corner, and the vector field was wavier than the baseline layout. This vortex flow in the x-direction is another flow characteristic attributed to the stacking of mixing unit in the cross flow direction. micromixer showed a meaningful enhancement of the mixing performance for Reynolds numbers less than about 20, as compared with other micromixers. The enhancement was caused by the two flow characteristics: the baffle impingement by a stream of high concentration and a swirl motion inside the mixing unit. The baffle impingement by a stream of high concentration was observed throughout all Reynolds numbers, while the swirl motion inside the mixing unit became significant only for Reynolds numbers larger than about five. Figure 13b shows that the required pressure load of the present micromixer was comparable with those of other micromixers; the mixing enhancement was achieved without an additional pressure load. This confirms that the mixing unit stacking in the cross flow direction is a practical approach to enhance the mixing performance, especially for Reynolds numbers less than about 20.

results, a rhombic angle of 90° provides the best mixing performance. The mixing performance was enhanced by the two pairs of counter rotating vortices in the cross section. The present

*Micromachines* **2021**, *12*, x FOR PEER REVIEW 14 of 18

**Figure 12.** Streamlines starting from the two inlets and concentration contours at the plane of the first baffle within each mixing unit for *Re* = 20. (**a**) Baseline. (**b**) Eight step stacking. **Figure 12.** Streamlines starting from the two inlets and concentration contours at the plane of the first baffle within each mixing unit for *Re* = 20. (**a**) Baseline. (**b**) Eight step stacking.

Figure 13 compares the simulated mixing performance in terms of the *DOM* and required pressure load with those from other passive micromixers. In general, the mixing performance of micromixers increased with the diffusion coefficient. The value of diffusion coefficient depends on the fluid used in actual applications. To neglect the effects of the

diffusion coefficient on the mixing performance, several passive micromixers based on the same size of the diffusion coefficient were selected for comparison; the diffusion coefficient was 1.2 × 10−9 m2/s. Tripathi et al. [23] numerically investigated the mixing performance of spiral micromixers and observed transverse flow due to the centrifugal effect. A pair of counter rotating vortices was found to enhance the mixing performance in the spiral micromixer. Raza et al. [29] proposed a SAR mixing unit combined with baffles in a curved channel and simulated its mixing performance. They showed that the mixing performance was better than their earlier version of SAR micro-mixer, and attributed the mixing enhancement to the combined effects of collision, Dean vortex, and secondary flows. Makhsuda et al. [52] proposed an optimal combination of the cross channel split and recombine (CC-SAR) and the mixing cell with baffles (MC-B). They showed that the mixing performance of the optimized micromixer was noticeably improved compared with a micromixer based on the CC-SAR alone. Hossain et al. [53] studied a micromixer with unbalanced three split rhombic sub-channels. According to their results, a rhombic angle of 90◦ provides the best mixing performance. The mixing performance was enhanced by the two pairs of counter rotating vortices in the cross section. The present micromixer showed a meaningful enhancement of the mixing performance for Reynolds numbers less than about 20, as compared with other micromixers. The enhancement was caused by the two flow characteristics: the baffle impingement by a stream of high concentration and a swirl motion inside the mixing unit. The baffle impingement by a stream of high concentration was observed throughout all Reynolds numbers, while the swirl motion inside the mixing unit became significant only for Reynolds numbers larger than about five. Figure 13b shows that the required pressure load of the present micromixer was comparable with those of other micromixers; the mixing enhancement was achieved without an additional pressure load. This confirms that the mixing unit stacking in the cross flow direction is a practical approach to enhance the mixing performance, especially for Reynolds numbers less than about 20. *Micromachines* **2021**, *12*, x FOR PEER REVIEW 15 of 18

**Figure 13.** Comparison of the simulated mixing performance with other micromixers. (**a**) *DOM* vs. *Re*. (**b**) Pressure load vs. *Re*. about five. **Figure 13.** Comparison of the simulated mixing performance with other micromixers. (**a**) *DOM* vs. *Re*. (**b**) Pressure load vs. *Re*.

A new passive micro-mixer with the mixing units stacked in the cross flow direction

The numerical solutions showed that the mixing unit stacking in the cross flow direction noticeably enhanced the mixing performance for Reynolds numbers less than about 20, as compared with other micromixers. The stacking layout would be different according to the Reynolds number. For Reynolds numbers less than about two, the two step stacking in the cross flow direction showed the best mixing performance. On the other hand, the eight step stacking showed the best mixing performance for Reynolds

The mixing enhancement due to the mixing unit stacking in the cross flow direction was caused by two different flow characteristics. One is the baffle impingement by a stream of high concentration. A stream of high concentration was observed to impinge on the baffle wall throughout all Reynolds numbers. It was observed even in the limited convection flow regime such as *Re* = 0.5. The other one is the swirl motion inside the mixing unit in the cross flow direction, and it was significant for Reynolds numbers larger than

verse directions. The mixing performance of the present micro-mixer was evaluated numerically by examining the *DOM* and *MEC*. The numerical simulation was carried out for

Reynolds numbers ranging from 0.5 to 50.

numbers larger than about two.

**6. Conclusions**

#### **6. Conclusions**

A new passive micro-mixer with the mixing units stacked in the cross flow direction was proposed and estimated. It consisted of eight mixing units, and each mixing unit contained four baffles; the baffles were staggered alternatively in the cross flow and transverse directions. The mixing performance of the present micro-mixer was evaluated numerically by examining the *DOM* and *MEC*. The numerical simulation was carried out for Reynolds numbers ranging from 0.5 to 50.

The numerical solutions showed that the mixing unit stacking in the cross flow direction noticeably enhanced the mixing performance for Reynolds numbers less than about 20, as compared with other micromixers. The stacking layout would be different according to the Reynolds number. For Reynolds numbers less than about two, the two step stacking in the cross flow direction showed the best mixing performance. On the other hand, the eight step stacking showed the best mixing performance for Reynolds numbers larger than about two.

The mixing enhancement due to the mixing unit stacking in the cross flow direction was caused by two different flow characteristics. One is the baffle impingement by a stream of high concentration. A stream of high concentration was observed to impinge on the baffle wall throughout all Reynolds numbers. It was observed even in the limited convection flow regime such as *Re* = 0.5. The other one is the swirl motion inside the mixing unit in the cross flow direction, and it was significant for Reynolds numbers larger than about five.

Since both halves of the present micro-mixer in the transverse direction are planar, they can be simply stacked to construct the present micromixer. The present micromixer performed better than other passive micro-mixers for low Reynolds numbers less than about 20. Therefore, it is expected to be useful as a necessary part of lab-on-a-chip devices and micro-total analysis systems.

**Author Contributions:** Simulation, M.J.; formal analysis, D.-J.K.; visualization, M.J.; writing original draft, D.-J.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by BOKUK© 2021 Research Grant.

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

#### **Nomenclature**

*σ* standard deviation

## **References**