Surface characteristics and damage distributions of diamond wire sawn wafers for silicon solar cells

Bhushan Sopori; Srinivas Devayajanam; Prakash Basnyat; Bhushan Sopori; Srinivas Devayajanam; Prakash Basnyat

doi:10.3934/matersci.2016.2.669

AIMS Materials Science

2016, Volume 3, Issue 2: 669-685. doi: 10.3934/matersci.2016.2.669

Previous Article Next Article

Research article Topical Sections

Surface characteristics and damage distributions of diamond wire sawn wafers for silicon solar cells

1.
National Renewable Energy Laboratory, Golden, CO 80401, USA
2.
New Jersey Institute of Technology, Newark, NJ 07102, USA

Received: 20 April 2016 Accepted: 08 June 2016 Published: 12 June 2016

This paper describes surface characteristics, in terms of its morphology, roughness and near-surface damage of Si wafers cut by diamond wire sawing (DWS) of Si ingots under different cutting conditions. Diamond wire sawn Si wafers exhibit nearly-periodic surface features of different spatial wavelengths, which correspond to kinematics of various movements during wafering, such as ingot feed, wire reciprocation, and wire snap. The surface damage occurs in the form of frozen-in dislocations, phase changes, and microcracks. The in-depth damage was determined by conventional methods such as TEM, SEM and angle-polishing/defect-etching. However, because these methods only provide local information, we have also applied a new technique that determines average damage depth over a large area. This technique uses sequential measurement of the minority carrier lifetime after etching thin layers from the surfaces. The lateral spatial damage variations, which seem to be mainly related to wire reciprocation process, were observed by photoluminescence and minority carrier lifetime mapping. Our results show a strong correlation of damage depth on the diamond grit size and wire usage.

Keywords:

Citation: Bhushan Sopori, Srinivas Devayajanam, Prakash Basnyat. Surface characteristics and damage distributions of diamond wire sawn wafers for silicon solar cells[J]. AIMS Materials Science, 2016, 3(2): 669-685. doi: 10.3934/matersci.2016.2.669

Related Papers:

[1]	Zhonghua Qiao, Xuguang Yang . A multiple-relaxation-time lattice Boltzmann method with Beam-Warming scheme for a coupled chemotaxis-fluid model. Electronic Research Archive, 2020, 28(3): 1207-1225. doi: 10.3934/era.2020066
[2]	Yunfei Yuan, Changchun Liu . Optimal control for the coupled chemotaxis-fluid models in two space dimensions. Electronic Research Archive, 2021, 29(6): 4269-4296. doi: 10.3934/era.2021085
[3]	Jiwei Jia, Young-Ju Lee, Yue Feng, Zichan Wang, Zhongshu Zhao . Hybridized weak Galerkin finite element methods for Brinkman equations. Electronic Research Archive, 2021, 29(3): 2489-2516. doi: 10.3934/era.2020126
[4]	Wenya Qi, Padmanabhan Seshaiyer, Junping Wang . A four-field mixed finite element method for Biot's consolidation problems. Electronic Research Archive, 2021, 29(3): 2517-2532. doi: 10.3934/era.2020127
[5]	Shanpu Gao, Yubo Li, Anping Wu, Hao Jiang, Feng Liu, Xinlong Feng . An intelligent optimization method for accelerating physical quantity reconstruction in computational fluid dynamics. Electronic Research Archive, 2025, 33(5): 2881-2924. doi: 10.3934/era.2025127
[6]	Nisachon Kumankat, Kanognudge Wuttanachamsri . Well-posedness of generalized Stokes-Brinkman equations modeling moving solid phases. Electronic Research Archive, 2023, 31(3): 1641-1661. doi: 10.3934/era.2023085
[7]	Ziqiang Wang, Chunyu Cen, Junying Cao . Topological optimization algorithm for mechanical-electrical coupling of periodic composite materials. Electronic Research Archive, 2023, 31(5): 2689-2707. doi: 10.3934/era.2023136
[8]	Jiangshan Wang, Lingxiong Meng, Hongen Jia . Numerical analysis of modular grad-div stability methods for the time-dependent Navier-Stokes/Darcy model. Electronic Research Archive, 2020, 28(3): 1191-1205. doi: 10.3934/era.2020065
[9]	Yusen Lin . Periodic measures of reaction-diffusion lattice systems driven by superlinear noise. Electronic Research Archive, 2022, 30(1): 35-51. doi: 10.3934/era.2022002
[10]	Ling Xue, Min Zhang, Kun Zhao, Xiaoming Zheng . Controlled dynamics of a chemotaxis model with logarithmic sensitivity by physical boundary conditions. Electronic Research Archive, 2022, 30(12): 4530-4552. doi: 10.3934/era.2022230

Abstract

1. Introduction

The chemotaxis refers to the collective motion of bacteria, cells or an organisms in response to an attractant gradient. There are numerous examples of chemotaxis in animal and insect ecology, biological and biomedical sciences [9]. Deep understanding of the behavior of chemotaxis phenomenon is of great significance to the biological and medical applications.

Various kinds of mathematical models have been developed by both experimentalists and theoreticians to describe the chemotactic phenomenon. In addition, to further describe and interpret the bioconvection phenomenon [14], a variety of coupled chemotaxis-fluid models have been proposed and studied in [5,8,24]. The coupled model is expressed as follows:

$\begin{align} \frac{{\partial n}}{{\partial t}} + \nabla \cdot (\mathbf{u}n) = {D_n}{\nabla ^2}n - \chi \nabla \cdot [r(c)n\nabla c], \end{align}$ $

(1)

$\begin{align} \frac{{\partial c}}{{\partial t}} + \nabla \cdot (\mathbf{u}c) = {D_c}{\nabla ^2}c - n\kappa r(c), \end{align}$ $

(2)

$\begin{align} \rho (\frac{{\partial \mathbf{u}}}{{\partial t}} + \mathbf{u} \cdot \nabla \mathbf{u}) = - \nabla p + \eta {\nabla ^2}\mathbf{u} + n{V_b}({\rho _b} - \rho )\mathbf{g}, \end{align}$ $

(3)

$\begin{align} \nabla \cdot \mathbf{u} = {\rm{0}}. \end{align}$ $

(4)

Here, $ \mathbf{u} $ is velocity, $ \rho $ is the fluid density, $ p $ is the hydrodynamic pressure, $ \eta $ is the coefficient of dynamic viscosity, $ \mathbf{g} $ is the gravity, $ n $ and $ c $ represent the concentrations of bacteria and oxygen, respectively. $ \kappa $ is the oxygen consumption rate, $ \chi $ is the chemotactic sensitivity, $ D_n $ and $ D_c $ are diffusion coefficients for bacteria and oxygen, respectively. $ V_b $ and $ \rho_b $ are the bacterial volume and density, respectively. $ r(c) $ is the dimensionless cut-off function.

It can be seen from the above system that the chemotaxis-fluid model involves highly nonlinearity coupled with chemotaxis, diffusion and convection, which brings great difficulties to get the analytical solution. Recently, several numerical methods have been developed to study the chemotaxis phenomenon, such as finite element, finite volume schemes [13,21], upwind-difference potentials method [10], fractional step algorithms [20,25], interior penalty/discontinuous Galerkin methods [11,12], and cell-overcrowding prevention models [1,7], etc. In the work by Chertock et al. [5], a coupled chemotaxis-fluid system has been numerically studied by a high-resolution vorticity-based hybrid finite-volume finite-difference scheme. In [22], the chemotaxis phenomenon in incompressible viscous fluid has been numerically studied by a compact finite difference scheme. The plume patterns were studied by solving a chemotaxis-diffusion-convection coupling system in [6], where the comparison of different systems of chemotaxis-diffusion-convection, Rayleigh-B$ \acute{e} $nard convection and the double diffusive convection were investigated. And it was indicated that although physical mechanisms are different, the convection patterns and the dimensionless differential systems are similar. In addition, a three-dimensional chemotaxis-fluid model was numerically investigated in [17]. The formation of falling bacterial plumes and its convergence to a steady state have been observed.

In this work, we will develop an efficient LBM to solve the coupled chemotaxis-fluid model. LBM originates from kinetic Boltzmann equation, which has been successfully applied in modeling complex fluid flows [18,4,2]. In addition, LBM also shows its excellent capability in solving the nonlinear systems [23,3,28,26]. For the chemotaxis problem, a simplified LBM was developed for a bacterial chemotaxis model by Hilpert [16]. While, Yan et al. [27] developed a MRT-LBM to investigate the traveling bacterial bands of chemotactic bacteria. However, due to the highly nonlinearity of the coupled chemotaxis-fluid model, we need to develop more accurate and stable LBM to solve this complex problem. In [19], a novel MRT-LBM was proposed to simulate multi-phase fluids with Peng-Robinson (P-R) equation of state. To enhance the stability of the MRT-LBM, B-W scheme is introduced in the discretization of the evolution equation. Numerical results showed that the P-R free energy model was solved precisely and spurious currents were eliminated effectively. Inspired by this, we will apply the MRT-LBM with B-W scheme to solve the coupled chemotaxis-fluid model.

The rest of this article is organized as follows. In Section Ⅱ, the MRT-LBM with B-W scheme for the chemotaxis-fluid model is proposed. In Section Ⅲ, several numerical experiments are carried out on different chemotaxis problems, including the classical K-S model, the coupled Navier-Stokes-Keller-Segel (NS-KS) model, and the bacterial biocovection problem. The paper ends with some conclusions in Section Ⅳ.

2. Numerical method

In this section, we will introduce the MRT-LBM with B-W scheme for the coupled chemotaxis-fluid model (1)-(4). To better investigate the coupled chemotaxis-fluid model, the following dimensionless form is taken:

$\begin{align*} \mathbf{x'} = \frac{\bf{x}}{L}, \quad t' = \frac{D_n}{L^2}t, \quad c' = \frac{c}{c_{air}}, \quad n' = \frac{n}{n_r}, \\ \mathbf{u'} = \frac{L}{D_n}\mathbf{u}, \quad p' = \frac{L^2}{\eta D_n}p, \quad \mathbf{g'} = \frac{\bf{g}}{g}, \end{align*}$ $

where $ n_r $ is the characteristic cell density, $ c_{air} $ is the oxygen concentration in the air, and $ L $ is a characteristic length. Through introducing the above scaled variables and dropping the prime notation, the following system can be obtained.

$\begin{align} \frac{{\partial n}}{{\partial t}} + \nabla \cdot (\mathbf{u}n) = {\nabla ^2}n - \alpha \nabla \cdot [r(c)n\nabla c], \end{align}$ $

(5)

$\begin{align} \frac{{\partial c}}{{\partial t}} + \nabla \cdot (\mathbf{u}c) = \delta {\nabla ^2}c - \beta r(c)n, \end{align}$ $

(6)

$\begin{align} \frac{{\partial \mathbf{u}}}{{\partial t}} + \mathbf{u} \cdot \nabla \mathbf{u} = - \nabla p + S_c{\nabla ^2}\mathbf{u} + \gamma S_c n \mathbf{z}. \end{align}$ $

(7)

$\begin{align} \nabla \cdot \mathbf{u} = {\rm{0}}, \end{align}$ $

(8)

where $ \mathbf{z} $ is the unit vector in the vertical direction. $ r(c) $ is regularized by

$ r(c) = \frac{1}{2}(1+ \frac{c-c^*}{\sqrt{(c-c*)^2}+\epsilon^2}), $

where $ c^* $ is a cut-off constant and $ \epsilon>0 $ is a constant close to zero.

The dimensionless parameters $ \alpha $, $ \beta $, $ \gamma $, $ \delta $, and the Schmidt number $ S_c $ are given below

$\begin{align*} \alpha = \frac{\chi c_{air}}{D_n}, \quad \beta = \frac{\kappa n_r L^2}{c_{air} D_n}, \quad \gamma = \frac{V_b n_r g(\rho_b-\rho)L^3}{\eta D_n}, \quad \delta = \frac{D_c}{D_n}, \quad S_c = \frac{\eta}{D_n \rho}. \end{align*}$ $

2.1. Lattice Boltzmann method for the incompressible Navier-Stokes equations

In this section, LBM with MRT collision operator is applied to solve the incompressible N-S equations (7)-(8). To capture the chemotaxis phenomenon more accurately, B-W scheme is introduced in the discretization of the evolution equation of MRT-LBM.

The discrete velocity Boltzmann equation with MRT collision operator is expressed as

$\begin{equation} \frac{{\partial {h_i}}}{{\partial t}} + {\hat{c}\mathbf{e}_i} \cdot \nabla {h_i} = - {\Lambda _{ij}}[{h_j} - h_j^{eq}] + {H_i}, \end{equation}$ $

(9)

where, $ h_i (\mathbf{x}, t) $ is the discrete distribution function of particle at site $ \mathbf{x} $ and time $ t $ moving with speed $ \hat{c} $ along the direction $ \mathbf{e}_i $. $ \{\bf{e}_i, i = 0, 1, \cdots, k-1\} $ is the set of discrete velocity directions, $ \hat{c} $ is the sound speed, $ {\Lambda _{ij}} $ is the collision matrix, $ h_i^{eq}(\mathbf{x}, t) $ is the equilibrium distribution function (EDF), and $ H_i $ is the force distribution function.

The discrete velocity Boltzmann equation (9) is solved by a time-splitting scheme, which is decomposed into two sub-processes, i.e., the collision process,

$\begin{equation} \frac{{\partial {h_i}}}{{\partial t}} = - {\Lambda _{ij}}[{h_j} - h_j^{eq}] + {H_i}, \end{equation}$ $

(10)

and the streaming process,

$\begin{equation} \frac{{\partial {h_i}}}{{\partial t}} + {\hat{c}\mathbf{e}_i} \cdot \nabla {h_i} = 0. \end{equation}$ $

(11)

In the MRT model, the collision subprocess is carried out in the moment space. We take the generally used D2Q9 (two-dimensional space with nine discrete velocities) model as an example. The distribution functions $ h_i $ in the moment space are defined as

$ \mathbf{m} = \mathbf{M} \cdot \mathbf{h} = (\rho, e, \varepsilon, j_x, q_x, j_y, q_y, p_{xx}, p_{xy})^T, $

where $ e $ and $ \varepsilon $ are related to total energy and the function of energy square; $ j_x $ and $ j_y $ are relevant to the momentum; $ q_x $ and $ q_y $ are related to the $ x $ and $ y $ components of the energy; $ p_{xx} $ and $ p_{xy} $ are the corresponding diagonal and off-diagonal components of the stress tensor, respectively. $ \mathbf{M} $ is the transformation matrix. $ \mathbf{h} = \{h_i, i = 0, 1, \cdots k-1\} $. For the D2Q9 model, $ \mathbf{M} $ is defined as

$ \mathbf{M} = \left( {

$\begin{array}{*{20}{c}} 1&1&1&1&1&1&1&1&1\\ { - 4}&{ - 1}&{ - 1}&{ - 1}&{ - 1}&2&2&2&2\\ 4&{ - 2}&{ - 2}&{ - 2}&{ - 2}&1&1&1&1\\ 0&1&0&{ - 1}&0&1&{ - 1}&{ - 1}&1\\ 0&{ - 2}&0&2&0&1&{ - 1}&{ - 1}&1\\ 0&0&1&0&{ - 1}&1&1&{ - 1}&{ - 1}\\ 0&0&{ - 2}&0&2&1&1&{ - 1}&{ - 1}\\ 0&1&{ - 1}&1&{ - 1}&0&0&0&0\\ 0&0&0&0&0&1&{ - 1}&1&{ - 1} \end{array}$ } \right). $

Through multiplying the transformation matrix $ \mathbf{M} $, the collision process can be rewritten onto the moment space as

$\begin{equation} \frac{{\partial {\bf{m}}}}{{\partial t}} = - \mathbf{\widetilde{S}}(\mathbf{m}-\mathbf{m}^{eq})+ \widehat{\bf{H}}, \end{equation}$ $

(12)

where $ \mathbf{\widetilde{S}} = \mathbf{M\Lambda M^{-1}} $ is a diagonal relaxation matrix, given by

$ \mathbf{\widetilde{S}} = \mathrm{diag}\{\widetilde{s}_0, \widetilde{s}_1, \widetilde{s}_2, \widetilde{s}_3, \widetilde{s}_4, \widetilde{s}_5, \widetilde{s}_6, \widetilde{s}_7, \widetilde{s}_8\}. $

The equilibrium moment $ \mathbf{m}^{eq} $ is defined as

$\begin{equation} {\mathbf{m}^{eq}} = \mathbf{M} \cdot {\mathbf{h}^{eq}} = \rho \left[ {\begin{array}{*{20}{c}} 1\\ { - 2 + 3{\mathbf{u}^2}}\\ {1 - 3{\mathbf{u}^2}}\\ u\\ { - u}\\ v\\ { - v}\\ {{u^2} - {v^2}}\\ {uv} \end{array}} \right]. \end{equation}$ $

(13)

$ \widehat{\bf{H}} = \mathbf{M}\mathbf{H} $ is the corresponding force moment, which is expressed as

$\begin{equation} \begin{array}{l} {\widehat{H}_0} = 0, \\ {\widehat{H}_1} = 6(1 - \frac{{{s_1}}}{2})\mathbf{u} \cdot \mathbf{F}_t, \\ {\widehat{H}_2} = - 6(1 - \frac{{{s_2}}}{2})\mathbf{u} \cdot \mathbf{F}_t, \\ {\widehat{H}_3} = {F_{tx}}, \\ {\widehat{H}_4} = - (1 - \frac{{{s_4}}}{2}){F_{tx}}, \\ {\widehat{H}_5} = {F_{ty}}, \\ {\widehat{H}_6} = - (1 - \frac{{{s_6}}}{2}){F_{ty}}, \\ {\widehat{H}_7} = 2(1 - \frac{{{s_7 }}}{2})(u{F_{tx}} - v{F_{ty}}), \\ {\widehat{H}_8} = (1 - \frac{{{s_8}}}{2})(u{F_{ty}} + v{F_{tx}}), \end{array} \end{equation}$ $

(14)

where $ \mathbf{H} = (H_0, ..., H_8)^T $. $ \mathbf{F}_t $ is the external force, which has the following form

$ \mathbf{F}_t = (F_{tx}, F_{ty}). $

In this simulation, $ \mathbf{F}_t $ should be expressed as

$\begin{equation} \mathbf{F}_t = \gamma S_c n \mathbf{z}. \end{equation}$ $

(15)

The first order explicit Euler scheme is used to discrete (12) as

$\begin{equation} \mathbf{m}^+ = \mathbf{m} - \mathbf{S}(\mathbf{m}-\mathbf{m}^{eq})+ \delta t\widehat{\bf{H}}, \end{equation}$ $

(16)

where $ \mathbf{m}^+ = \mathbf{M} \mathbf{h}^+ $ is the post-collision moments, and $ \mathbf{h}^+ = (h^+_0, ..., h^+_8)^T $ is the post-collision distribution function. The non-dimensional relaxation matrix is defined as $ \mathbf{S} = \delta t \mathbf{\widetilde{S}} $, which is given by

$ \mathbf{S} = \mathrm{diag}\{s_0, s_1, s_2, s_3, s_4, s_5, s_6, s_7, s_8\}. $

In simulations, $ s_0 = s_3 = s_5 $, $ s_4 = s_6 $ and $ s_7 = s_8 $.

Next, through using the second-order B-W scheme, Eq. (11) can be solved on a regular lattice with spacing $ \delta x $,

$\begin{equation} \begin{split} {h_i}(\mathbf{x}, t + \delta t) = & h_i^ + (\mathbf{x}, t) - \frac{A}{2}(3 h_i^ + (\mathbf{x}, t) - 4 h_i^ + (\mathbf{x} - {\bf{e}_i}\delta x, t) + h_i^ + (\mathbf{x} - 2{\bf{e}_i}\delta x, t))\\ & + \frac{A^2}{2}(h_i^ + (\mathbf{x}, t) - 2 h_i^ + (\mathbf{x} - {\bf{e}_i}\delta x, t) + h_i^ + (\mathbf{x} - 2{\bf{e}_i}\delta x, t)), \end{split} \end{equation}$ $

(17)

where the time step $ \delta t $ is related to the CFL number $ A $ as $ \delta t = A \delta x /c $.

The macroscopic quantities, $ \rho $ and $ \mathbf{u} $ are calculated as

$\begin{equation} \rho = \sum\limits_i {{h_i}} , \quad\rho \mathbf{u} = \sum\limits_i {{\bf{c}_i}{h_i}} + \frac{{\delta t}}{2}\mathbf{F}_t. \end{equation}$ $

(18)

2.2. Lattice Boltzmann method for the transport equation of concentrations of bacteria and oxygen

The evolution equations of the LBM for (5) and (6) are given by

$\begin{align} \frac{{\partial {f_i}}}{{\partial t}} + {\bf{\hat{c}}_i} \cdot \nabla {f_i} = - \frac{1}{{{\tau _n}}}({f_i} - f_i^{eq}) + {S_i}, \end{align}$ $

(19)

$\begin{align} \frac{{\partial {g_i}}}{{\partial t}} + {\bf{\hat{c}}_i} \cdot \nabla {g_i} = - \frac{1}{{{\tau _c}}}({g_i} - g_i^{eq}) + {\hat{S}_i} + \frac{{\delta t}}{2}\partial t{\hat{S}_i}, \end{align}$ $

(20)

where $ \tau_n $ and $ \tau_c $ are the dimensionless relaxation factors, $ \mathbf{\hat{c}}_i = \hat{c}\mathbf{e}_i $, $ f_i ({\bf{x}}, t) $ and $ g_i ({\bf{x}}, t) $ are the distribution functions of particle moving with velocity $ \mathbf{\hat{c}}_i $ at position $ \mathbf{x} $ and time $ t $. $ S_i ({\bf{x}}, t) $ and $ \hat{S}_i({\bf{x}}, t) $ are source distribution functions, which are defined as

$\begin{align*} {S_i} = {\omega _i}{\bf{\hat{c}}_i} \cdot B r(c)n\nabla c, \quad \hat{S}_i = {\omega _i}( - \beta r(c)n). \end{align*}$ $

Here, $ B $ is a parameter, which satisfies $ \alpha = B c_s^2{\tau _n}\delta t $. $ f_i^{eq} (\mathbf{x}, t) $ and $ g_i^{eq} (\mathbf{x}, t) $ are local equilibrium distribution functions, and read as

$\begin{align*} f_i^{eq}(\mathbf{x}, t) = {\omega _i}n(1 + \frac{{{\bf{\hat{c}}_i} \cdot \mathbf{u}}}{{c_s^2}}), \quad g_i^{eq}(\mathbf{x}, t) = {\omega _i}c(1 + \frac{{{\bf{\hat{c}}_i} \cdot \mathbf{u}}}{{c_s^2}}). \end{align*}$ $

(19) and (20) are decomposed into two sub-processes: the collision process

$\begin{align} \frac{{\partial {f_i}}}{{\partial t}} = - \frac{1}{{{\tau _n}}}({f_i} - f_i^{eq}) + {S_i}, \end{align}$ $

(21)

$\begin{align} \frac{{\partial {g_i}}}{{\partial t}} = - \frac{1}{{{\tau _c}}}({g_i} - g_i^{eq}) + {\hat{S}_i} + \frac{{\delta t}}{2}\partial t{\hat{S}_i}, \end{align}$ $

(22)

and the streaming process,

$\begin{align} \frac{{\partial {f_i}}}{{\partial t}} + {\mathbf{\hat{c}}_i} \cdot \nabla f_i = 0, \end{align}$ $

(23)

$\begin{align} \frac{{\partial {g_i}}}{{\partial t}} + {\mathbf{\hat{c}}_i} \cdot \nabla {g_i} = 0 . \end{align}$ $

(24)

The first order explicit Euler scheme is used to solve the collision process (21)-(22), which leads to

$\begin{align*} f_i^ + = {f_i} - \frac{{\delta t}}{{{\tau _n}}}({f_i} - f_i^{eq}) + \delta t{S_i}, \\ g_i^ + = {g_i} - \frac{{\delta t}}{{{\tau _c}}}({g_i} - g_i^{eq}) + \delta t{\hat{S}_i} + \frac{{\delta {t^2}}}{2}\partial t{\hat{S}_i}. \end{align*}$ $

B-W scheme is used to solve (23) and (24),

$\begin{equation} \begin{split} {f_i}(\mathbf{x}, t + \delta t) = & f_i^ + (\mathbf{x}, t) - \frac{A}{2}(3 f_i^ + (\mathbf{x}, t) - 4 f_i^ + (\mathbf{x} - {\mathbf{e}_i}\delta x, t) + f_i^ + (\mathbf{x} - 2{\mathbf{e}_i}\delta x, t))\\ & + \frac{A^2}{2}(f_i^ + (\mathbf{x}, t) - 2 f_i^ + (\mathbf{x} - {\mathbf{e}_i}\delta x, t) + f_i^ + (\mathbf{x} - 2{\mathbf{e}_i}\delta x, t)), \end{split} \end{equation}$ $

(25)

$\begin{equation} \begin{split} {g_i}(\mathbf{x}, t + \delta t) = & g_i^ + (\mathbf{x}, t) - \frac{A}{2}(3 g_i^ + (\mathbf{x}, t) - 4 g_i^ + (\mathbf{x} - {\mathbf{e}_i}\delta x, t) + g_i^ + (\mathbf{x} - 2{\mathbf{e}_i}\delta x, t))\\ & + \frac{A^2}{2}(g_i^ + (\mathbf{x}, t) - 2 g_i^ + (\mathbf{x} - {\mathbf{e}_i}\delta x, t) + g_i^ + (\mathbf{x} - 2{\mathbf{e}_i}\delta x, t)). \end{split} \end{equation}$ $

(26)

The concentrations of the bacteria $ n(\mathbf{x}, t) $ and oxygen $ c(\mathbf{x}, t) $ are computed by

$\begin{align*} n (\mathbf{x}, t) = \sum\limits_i {f_i (\mathbf{x}, t)}, \quad c (\mathbf{x}, t) = \sum\limits_i {g_i (\mathbf{x}, t)}. \end{align*}$ $

3. Numerical experiments

To validate that the stability is improved by the present LBM with B-W scheme, and also to verify the capability of the proposed method, a series of chemotaxis problems are numerically investigated. In the following numerical experiments, the general bounce-back scheme in reference [29] is used to treat the macroscopic boundary conditions in the LBM with B-W scheme.

3.1. Validation of the LBM with B-W scheme

In this section, the following classical K-S chemotaxis problem is considered.

$\begin{align} {{{\partial n} \over {\partial t}} = \mu \nabla \cdot n - \chi \nabla \cdot (n \nabla c), } \end{align}$ $

(27)

$\begin{align} {{{\partial c} \over {\partial t}} = \nabla ^2 c + n - c.} \end{align}$ $

(28)

(27)-(28) can be regarded as a special case of the generalized K-S model. $ \mu $ and $ \chi $ represent the diffusion constants. Without loss of generality, the values of these parameters are set as $ \mu = \chi = 1 $. The computational domain is $ \Omega = [-0.5, 0.5]\times[-0.5, 0.5] $ and the initial conditions are

$ n(x, y, 0) = 1000e^{ - 100(x^2 + y^2 )}, \quad c(x, y, 0) = 500e^{ - 50(x^2 +y^2 )}. $

The boundary conditions are

$ {{\partial n} \over {\partial \mathbf{v}}} = {{\partial c} \over {\partial \mathbf{v}}} = 0, (x, y) \in \partial \Omega, $

where $ \mathbf{v} $ is unit normal vector and $ \partial \Omega $ represents the domain boundary.

Under above initial and boundary conditions, the solution of (27)-(28) will blow up at the origin in finite time, which brings a great challenge in the numerical simulation. To illustrate that the stability is improved by B-W scheme, the comparison between standard LBM and the present LBM with B-W scheme is carried out. In the numerical experiment, the D2Q9 lattice model is used with a $ 256\times256 $ uniform grid in space and time step $ \delta t = 5.0 \times 10^{-7} $. Fig. 1 shows numerical solutions of bacteria concentrations $ n(\mathbf{x}, t) $ at different time points by the standard LBM, i.e., $ A = 0 $. It can be observed that numerical solutions appear negative values as time evolves. Then, the LBM with B-W scheme is applied to solve the above initial-value problem. In the simulation, the same grid is used with the CFL number $ A = 0.5 $, which gives $ \delta t = 2.5 \times 10^{-7} $. Numerical solutions can be found in Fig. 2, which show that the positivity of bacteria concentrations $ n(\mathbf{x}, t) $ is well preserved. In addition, one-dimensional profiles of $ n(\mathbf{x}, t) $ along $ y = 0 $ at $ t = 1.15\times10^{-4} $ are depicted in Fig. 3. Negative values can be observed clearly near the origin in Fig. 3(a), while not in Fig. 3(b), which verifies that the numerical stability is improved by the proposed LBM with B-W scheme.

Figure 1. Numerical solutions of bacteria concentrations $ n(\mathbf{x}, t) $ of (27)-(28) by the standard LBM (A = 0): (a) $ t = 5.0\times10^{-6} $, (b) $ t = 1.0\times10^{-5} $, (c) $ t = 4.5\times10^{-5} $, (d) $ t = 1.15\times10^{-4} $.

DownLoad: Full-Size Img PowerPoint

Figure 2. Numerical solutions of bacteria concentrations $ n(\mathbf{x}, t) $ of (27)-(28) by the present LBM with B-W scheme (A = 0.5), (a) $ t = 5.0\times10^{-6} $, (b) $ t = 1.0\times10^{-5} $, (c) $ t = 4.5\times10^{-5} $, (d) $ t = 1.15\times10^{-4} $.

DownLoad: Full-Size Img PowerPoint

Figure 3. 1D profiles of $ n(\mathbf{x}, t) $ along y = 0 at $ t = 1.15\times10^{-4} $: (a) by the standard LBM (A = 0); (b) by the present LBM with B-W scheme (A = 0.5).

DownLoad: Full-Size Img PowerPoint

3.2. Numerical study of the coupled NS-KS equations

In this section, the incompressible N-S equations coupled with K-S equation is solved for the convergence test of the proposed LBM with B-W scheme. The mathematical model has the following form,

$\begin{align*} \frac{{\partial n}}{{\partial t}} + \nabla \cdot (\mathbf{u}n) = {\nabla ^2}n - \nabla \cdot [\chi n\nabla c], \\ \frac{{\partial c}}{{\partial t}} + \nabla \cdot (\mathbf{u}c) = {\nabla ^2}c + n f(c) - c\kappa, \end{align*}$ $

$\begin{align} \frac{{\partial \mathbf{u}}}{{\partial t}} + \mathbf{u} \cdot \nabla \mathbf{u} = - \nabla p + \mu {\nabla ^2}\mathbf{u}, \end{align}$ $

(29)

$\begin{align} \nabla \cdot \mathbf{u} = {\rm{0}}. \end{align}$ $

(30)

In the numerical simulation, $ f = \kappa = 1 $. The boundary conditions are

$\begin{align*} {{\partial n} \over {\partial \mathbf{v}}} = {{\partial c} \over {\partial \mathbf{v}}} = {{\partial p} \over {\partial \mathbf{v}}} = 0, \quad \mathbf{u} = \mathbf{u}_{exact}, \quad (x, y) \in \partial \Omega. \end{align*}$ $

The computational domain is $ \Omega = [0, 1]\times[0, 1] $. The above problem has the following exact solutions [22],

$\begin{align*} &u_{exact} = -\cos(\pi x)\sin(\pi y) e^{-2\pi^2 t}, \quad v_{exact} = -\sin(\pi x)\cos(\pi y) e^{-2\pi^2 t}, \\ &p = c_1-0.25(\cos(2\pi x) + \cos(2\pi y)) e^{-4\pi^2 t}, \\ &n_{exact} = \cos(\pi x)\cos(\pi y) e^{-\pi^2t}, \quad c_{exact} = \cos(\pi x)\cos(\pi y) e^{-2\pi^2 t}. \end{align*}$ $

To measure the accuracy, the following relative error is applied:

$\begin{equation*} E = \frac{\Sigma |\phi(\mathbf{x}, t) - \phi^*(\mathbf{x}, t)|}{\Sigma |\phi^*(\mathbf{x}, t)|}, \end{equation*}$ $

where $ \phi^* $ and $ \phi $ represent the analytical and numerical solutions, respectively, and the summation is taken over all grid points.

It can be seen from Table 1 that numerical solutions of $ n, c, \mathbf{u} $ by the proposed LBM with B-W scheme has second order convergence.

Table 1. Convergence test with different lattice spacings, $ \delta t = 1.0\times10^{-5} $.

$ Meshes $	$ E_n $	order	$ E_c $	order	$ E_\mathbf{u} $	order
$ 64\times 64 $	$ 1.3934 \times 10^{-4} $	–	$ 1.6327 \times 10^{-4} $	–	$ 2.0029 \times 10^{-4} $	–
$ 128\times 128 $	$ 3.6892 \times 10^{-5} $	1.9172	$ 4.3355 \times10^{-5} $	1.9130	$ 5.3150 \times 10^{-5} $	1.9139
$ 256\times 256 $	$ 9.4955 \times 10^{-6} $	1.9580	$ 1.1174\times 10^{-5} $	1.9560	$ 1.3699 \times 10^{-5} $	1.9560
$ 512\times 512 $	$ 2.4090 \times 10^{-6} $	1.9788	$ 2.8380\times 10^{-6} $	1.9772	$ 3.4775 \times 10^{-6} $	1.9779

| Show Table

DownLoad: CSV

3.3. Numerical study of the coupled chemotaxis-fluid system

The chemotaxis response of bacterial suspensions can be described by the chemotaxis-fluid model (5)-(8). The balance between chemotaxis, diffusion, and convection of bacteria leads to the formation and stability of plumes. There still needs to a deep understanding of the particular impact of each mechanism.

In the simulations, the parameters of (5)-(8) are set as $ \gamma = 418 $, $ S_c = 7700 $, $ c^* = 0.3 $, the grid spacing and time step are $ h = 1/128 $ and $ \delta t = 2\times10^{-5} $, respectively, and $ \epsilon = h $. The initial conditions are defined as:

$\begin{align*} {n^0}(x, y) & = \left\{ {\begin{array}{*{20}{l}} {1, \quad \quad{\rm{ if }} \quad y > 0.501 - 0.01(\sin (x - 0.5)\pi ))}, \\ {0.5, \quad {\rm{ otherwise, }}} \end{array}} \right.\\ c^0(x, y) & = 1, \quad \mathbf{u}^0(x, y) = 0, \end{align*}$ $

At the top boundary $ \Omega_{top} $, the boundary conditions are given as

$\begin{equation} \chi n r(c) c_y - D_n n_y = 0, c = 1.0, v = 0, u_y = 0, \forall (x, y) \in \partial \Omega_{top}, \end{equation}$ $

(31)

where $ \mathbf{u} = (u, v) $. At the bottom of the domain $ \partial \Omega_{bot} $, the boundary conditions are:

$\begin{equation} n_y = c_y = 0, u = v = 0, \forall (x, y) \in \partial \Omega_{bot}. \end{equation}$ $

(32)

The mixed boundary condition in Eq. (31) can be regrouped as $ \alpha c_y = (\ln n)_y $, where $ \alpha = \chi r(c)/D_n $. By integrating $ \alpha c_y = (\ln n)_y $ with respect to y, the mixed boundary condition is transformed into Dirichlet boundary condition, then the general bounce-back scheme in reference [29] is used.

First of all, as the numerical experiments in [17], the case of $ \alpha = 5 $, $ \beta = 5 $ and $ \delta = 0.25 $ is considered. Numerical results of the vertical profiles of $ c $ and cell densities $ n $ at t = 0.22 are shown in Fig. 4(a), which agree well with the available data in [17]. As a result of the cut-off of the chemotactic convection for oxygen levels below $ c\leq c^* $, an increase of cells towards the bottom can be found from the vertical profile of the cell density. The streamlines in Fig. 4(b) shows that the bioconvection phenomenon is well captured. Next, the effect of $ \beta $ is investigated for $ \alpha = 5 $ and $ \delta = 1 $. Numerical results for $ \beta = 7.23 $, 10, 20, 40 are shown in Fig. 5. When the value of $ \beta $ increases, the oxygen consumption increases relative to the oxygen diffusion. It can be seen from Fig. 5(a) that the oxygen density decreases with the increases of the value of $ \beta $, while the cell up-swimming increases (see Fig. 5(b)). The effect of $ \alpha $ is also investigated for the fixed value of $ \beta = 10 $ and $ \delta = 1 $. The cases of $ \alpha = 1 $, 2, 4, 5.952 are considered. Numerical results are shown in Fig. 6, which indicate that the cell density near the surface increases when the value of $ \alpha $ increases (see Fig. 6(b)), and less overall oxygen consumption occurs in these regions (see Fig. 6(a)). Similar results as in Fig. 5 and Fig. 6 can be found in [17].

Figure 4. (a) Vertical profiles of the oxygen $ c $ and cell densities $ n $ and (b) streamlines at t = 0.22.

DownLoad: Full-Size Img PowerPoint

Figure 5. Vertical profiles of (a) oxygen $ c $ and (b) cell densities $ n $ at t = 0.22 with $ \alpha = 5 $, $ \delta = 1 $ and different $ \beta $.

DownLoad: Full-Size Img PowerPoint

Figure 6. Vertical profiles of (a) oxygen $ c $ and (b) cell densities $ n $ at t = 0.22 with $ \beta = 10 $, $ \delta = 1 $ and different $ \alpha $.

DownLoad: Full-Size Img PowerPoint

3.4. Numerical study of bacterial bioconvection

In this subsection, the system (5)-(8) on the rectangular domain is numerically studied. In the following numerical simulations, the coefficients are set as $ \alpha = 10 $, $ \beta = 10 $, $ \gamma = 1000 $, $ \delta = 5 $, $ S_c = 500 $, and the cut-off $ c^* = 0.3 $. The computational domain is $ \Omega = [0, 6]\times[0, 1] $. The mesh grid and time step are $ h = 1/256 $ and $ \delta t = 10^{-5} $, respectively, and we use $ \epsilon = h $. The boundary conditions (31) and (32) are applied.

3.4.1. The case of randomly perturbed homogeneous initial data

As in [17], the following initial data is considered:

$\begin{align*} n^0(x, y) & = 0.8+0.2 {\rm{rand()}}, c^0(x, y) = 1, \\ \mathbf{u}^0(x, y) & = 0, \end{align*}$ $

where rand() is a random number uniformly distributed in the interval [0, 1].

The evolution of the cell density $ n $ at different time is depicted in Fig. 7. It can be observed that the cell concentration increases with time as a result of the cells accumulating at the top of the chamber, while the distinct phenomenon of the cell concentration can be found at the bottom of the chamber. As a consequence, several instabilities at a high-concentration layer is triggered below the fluid-air surface at $ t = 0.32 $. At $ t = 0.37 $, two falling plumes is formulated due to the instabilities and all plumes hit the bottom of the chamber at $ t = 0.42 $. At the bottom of the plumes, the cells move onto the edge of the plumes, and are advected by the flow field near the fluid-air surface. Thus the smaller plumes and thicker cell-rich layer at $ t = 0.47 $ than those at $ t = 0.42 $ can be observed. The velocity components $ u(x, y, t) $ and $ v(x, y, t) $ are depicted as arrows in Fig. 7. It can be seen that the fluid velocity is downward in regions of large concentrations. The above phenomenon was also observed in [5,17,6].

Figure 7. Bioconvection patterns of the cell density $ n(x, y) $ at (a) $ t = 0.32 $, (b) $ t = 0.37 $, (c) $ t = 0.42 $, and (d) $ t = 0.47 $.

DownLoad: Full-Size Img PowerPoint

3.4.2. Effect of increased value of $ \delta $

In this subsection, the effect of the parameter $ \delta $ on evolution of the cell density is investigated. In numerical simulations, the following initial data is used.

$\begin{align*} {n^0}(x, y) & = \left\{ {\begin{array}{*{20}{l}} {1, \quad \quad {\rm{ if }} \quad y > 0.499 - 0.05(\sin (x - 0.5)\pi ))}\\ {0.5, \quad {\rm{ otherwise, }}} \end{array}} \right.\\ c^0(x, y) & = 1, \quad \mathbf{u}^0(x, y) = \mathbf{0}. \end{align*}$ $

The computational domain is $ \Omega = [0, 4]\times[0, 1] $ and the values of $ \delta $ are set as: $ \delta = 5 $, 25, 50, while the value of $ D_n $ is fixed. The other parameters are $ \alpha = 10 $, $ \beta = 10 $, $ \gamma = 1000 $, and $ S_c = 500 $.

Numerical results of $ \delta = 5, $ 25, and 50 are shown in Figs. 8-10, which are consistent with the results in [17,15]. For $ \delta = 5 $ (see Fig. 8), profiles of plumes are similar to those in Fig. 7. However, for $ \delta = 25 $ and 50 (see Figs. 9 and 10, respectively), profiles of plumes are significantly different, where we could observe that the two instabilities are approaching each other and merge into one instability.

Figure 8. For $ \delta = 5 $, cell concentration $ n(x, y) $ at (a) $ t = 0.27 $, (b) $ t = 0.34 $, and (c) $ t = 0.39 $.

DownLoad: Full-Size Img PowerPoint

Figure 9. For $ \delta = 25 $, cell concentration $ n(x, y) $ at (a) $ t = 2.9 $, (b) $ t = 3.4 $, and (c) $ t = 3.9 $.

DownLoad: Full-Size Img PowerPoint

Figure 10. For $ \delta = 50 $, cell concentration $ n(x, y) $ at (a) $ t = 1.9 $, (b) $ t = 2.4 $, and (c) $ t = 2.9 $.

DownLoad: Full-Size Img PowerPoint

4. Conclusions

An accurate and stable MRT-LBM with B-W scheme for a chemotaxis-fluid model is developed. Through introducing B-W scheme in the evolution process of the LBM, the numerical stability is improved. The numerical study of the classical K-S model shows that the proposed LBM with B-W scheme could preserve the positivity of bacteria concentrations. Then, the second order accuracy is tested by simulating the coupled NS-KS model. The nonlinear dynamics of the chemotaxis-fluid model was investigated numerically by the proposed MRT-LBM with B-W scheme. Our numerical results agree well with those in the literature under different settings. In the numerical simulation of bioconvection phenomenon, the evolution of the instability of falling plumes and the convergence towards numerically stable stationary plumes are observed.

Appendix A. From MRT-LBE to the Hydrodynamic equations: Multi-scale Chapman-Enskog expansion

Rewritten the evolution equation (17) up to $ O(\delta x^3) $, one can obtain the following equation

$\begin{equation} {h_i}(\mathbf{x}, t + \delta t) = h_i^ + (\mathbf{x}, t) - \delta t c{\mathbf{e}_i} \cdot \nabla h_i^ + (\mathbf{x}, t) + \frac{1}{2}\delta {t^2}{(c{\mathbf{e}_i} \cdot \nabla )^2}h_i^ + (\mathbf{x}, t) + O(\delta {x^3}). \end{equation}$ $

(33)

Multiplying Eq. (16) by inverse of the transformation matrix $ M $, and substituting it into the above equation and expanding the variables around $ (\mathbf{x}, t) $ up to $ O(\delta x^2) $ and $ O(\delta t^2) $, we obtain

$\begin{equation} {\partial _t}{h_i} + c{\mathbf{e}_i} \cdot \nabla {h_i} = {\Omega _i} - \frac{{\delta t}}{2}[\partial _t^2{h_i} - {(c{\mathbf{e}_i} \cdot \nabla )^2}{h_i} + 2c{\mathbf{e}_i} \cdot \nabla {\Omega _i}] + O(\delta {x^2} + \delta {t^2}), \end{equation}$ $

(34)

where $ \Omega_i = - {\Lambda _{ij}}({h_j} - h_j^{eq}) + {H_i} $.

$ \partial _t^2{h_i} $ can be expressed as

$\begin{equation} \partial _t^2{h_i} = {(c{\mathbf{e}_i} \cdot \nabla )^2}{h_i} + {\partial _t}{\Omega_i} - c{\mathbf{e}_i} \cdot \nabla {\Omega _i} + O(\delta t). \end{equation}$ $

(35)

Thus, Eq. (34) can be rewritten as

$\begin{equation} {D_i}{h_i} = (1 - \frac{{\delta t}}{2}{D_i}){\Omega _i} + O(\delta {x^2} + \delta {t^2}), \end{equation}$ $

(36)

where $ D_i = {\partial _t} + c{\mathbf{e}_i} \cdot \nabla $.

In addition, from Eq. (34), we can see that $ {\Omega _i} = {D_i}{h_i} + O(\delta t) $. Thus, Eq. (36) is also equivalent to

$\begin{equation} {D_i}{h_i} + \frac{{\delta t}}{2}{D_i^2} h_i = - {\Lambda _{ij}}[{h_j} - h_j^{eq}] + {H_i} + O(\delta {x^2} + \delta {t^2}). \end{equation}$ $

(37)

Then we introduce the following expansions:

$\begin{equation} \begin{split} h_i & = h_i^{(0)} + \varepsilon h_i^{(1)} + \varepsilon ^2 h_i^{(2)} + \cdots, \\ {\partial \over {\partial t}}& = \varepsilon {\partial \over {\partial t_1}} + \varepsilon ^2 {\partial \over {\partial t_2}}, \nabla = \varepsilon \nabla_1, H_i = \varepsilon H_i^{(1)}, \end{split} \end{equation}$ $

(38)

where $ \varepsilon $ is a small parameter. With the expansions, Eq. (37) can be rewritten in consecutive orders of $ \varepsilon $ as

$\begin{equation} O({\varepsilon ^0}):{\rm{ }}h_i^{\left( 0 \right)} = h_i^{\left( {eq} \right)}, \end{equation}$ $

(39a)

$\begin{equation} O({\varepsilon ^1}):{\rm{ }}{D_{1i}}h_i^{\left( {eq} \right)} = - {\Lambda _{ij}}h_j^{(1)}+ H_i^{(1)}, \end{equation}$ $

(39b)

$\begin{equation} O({\varepsilon ^1}):{\rm{ }}{\partial _{{t_2}}}h_i^{(0)} + {D_{1i}}[({I_{ij}} - \frac{{{\Lambda _{ij}}}}{2})h_j^{(1)}] = - {\Lambda _{ij}}h_j^{(2)} - \frac{{\delta t}}{2}{D_{1i}}H_i^{(1)}, \end{equation}$ $

(39c)

where $ {D_{1i}} = {\partial _{{t_1}}} + {{\mathbf{c}}_i} \cdot {\nabla _1} $.

Multiplying the transformation Matrix $ \mathbf{M} $ on both side of Eq. (39), we can obtain the following moment equations:

$\begin{equation} O({\varepsilon ^0}):{\rm{ }}\mathbf{m}^{\left( 0 \right)} = \mathbf{m}^{\left( {eq} \right)}, \end{equation}$ $

(40a)

$\begin{equation} O({\varepsilon ^1}):{\rm{ }}\mathbf{D}_1{\mathbf{m}^{(0)}} = - \mathbf{\widetilde{S}}{\mathbf{m}^{(1)}} + {\mathbf{\widehat{H}}^{(1)}}, \end{equation}$ $

(40b)

$\begin{equation} O({\varepsilon ^2}):{\rm{ }}{\partial_{{t_2}}}{\mathbf{m}^{(0)}} + {\mathbf{D}_1}(\mathbf{I} - \frac{\mathbf{S}}{2}){\mathbf{m}^{(1)}} + \frac{{\delta t}}{2}{\mathbf{D}_1}{\mathbf{\widehat{H}}^{(1)}} = - \mathbf{\widetilde{S}}{\mathbf{m}^{(2)}}, \end{equation}$ $

(40c)

where $ \mathbf{D}_1 = \partial_{t_1}\mathbf{I}+\mathbf{C}_\alpha \partial_{1\alpha} $, $ \mathbf{C}_\alpha $ is the discrete velocity matrix.

In addition, from Eqs. (18) and (40a), we derive

$\begin{equation} \begin{split} \rho^{(1)} = 0, j_x^{(1)} = -\frac{\delta t}{2}F_{tx}^{(1)}, j_y^{(1)} = -\frac{\delta t}{2}F_{ty}^{(1)}, \rho^{(k)} = j_x^{(k)} = j_y^{(k)} = 0, k > 1. \end{split} \end{equation}$ $

(41)

On the $ t_1 $ time scale, Eq. (40b) can be rewritten as follows:

$\begin{align} &{\partial _{{t_1}}}\left[ {\begin{array}{*{20}{c}} \rho \\ {\rho ( - 2 + 3{\mathbf{u}^2})}\\ {\rho (1 - 3{\mathbf{u}^2})}\\ {\rho {u}}\\ { - \rho {u}}\\ {\rho {v}}\\ { - \rho {v}}\\ {\rho (u^2 - v^2)}\\ {\rho {u}{v}} \end{array}} \right] + {\partial _{1x}}\left[ {\begin{array}{*{20}{c}} {\rho {u}}\\ 0\\ { - \rho {u}}\\ {c_s^2\rho + \rho u^2}\\ {\rho {B_x}/3}\\ {\rho {u}{v}}\\ {\rho {u}{v}}\\ {2\rho {u}/3}\\ {\rho {v}/3} \end{array}} \right] + {\partial _{1y}}\left[ {\begin{array}{*{20}{c}} {\rho {v}}\\ 0\\ { - \rho {v}}\\ {\rho {u}{v}}\\ {\rho {u}{v}}\\ {c_s^2\rho + \rho v^2}\\ {\rho {B_y}/3}\\ { - 2\rho {v}/3}\\ {\rho {u}/3} \end{array}} \right] \\ = &\left[ {\begin{array}{*{20}{c}} 0\\ { - {\widetilde{s}_1}{e^{(1)}}}\\ { - {\widetilde{s}_2}{\varepsilon ^{(1)}}}\\ 0\\ { - {\widetilde{s}_4}q_x^{(1)}}\\ 0\\ { - {\widetilde{s}_6}q_y^{(1)}}\\ { - {\widetilde{s}_7}p_{xx}^{(1)}}\\ { - {\widetilde{s}_8}p_{xy}^{(1)}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} 0\\ {6(1 - {s_1}/2)\mathbf{u} \cdot {\mathbf{F}_t^{(1)}}}\\ { - 6(1 - {s_2}/2)\mathbf{u} \cdot {\mathbf{F}_t^{(1)}}}\\ {{F_{tx}^{(1)}}}\\ { - (1 - {s_4}/2){F_{tx}^{(1)}}}\\ {{F_{ty}^{(1)}}}\\ { - (1 - {s_6}/2){F_{ty}^{(1)}}}\\ {2(1 - {s_7}/2)(u{F_{tx}^{(1)}} - v{F_{ty}^{(1)}})}\\ {(1 - {s_8}/2)(u{F_{ty}^{(1)}} + v{F_{tx}^{(1)}})} \end{array}} \right] \end{align}$ $

(42)

where $ B_x = -1+6v^2+3\mathbf{u}^2 $, $ B_y = -1+6u^2+3\mathbf{u}^2 $.

Similarly, From Eq. (40b), the scale equations of conserved quantities $ \rho $, $ j_x $ and $ j_y $ on the $ t_2 $ time scale can be rewritten as

$\begin{equation} {\partial _{{t_2}}}\rho = 0. \end{equation}$ $

(43)

$\begin{equation} \begin{split} {\partial _{t2}}(\rho u) + \frac{1}{6}(1 - \frac{{{s_1}}}{2}){\partial _{1x}}{e^{(1)}} + (1 - \frac{{{s_7}}}{2})(\frac{1}{2}{\partial _{1x}}p_{xx}^{(1)} + {\partial _{1y}}p_{xy}^{(1)}) \\ + \frac{{\delta t}}{2}(1 - \frac{{{s_1}}}{2}){\partial _{1x}}(\mathbf{u} \cdot \mathbf{F}_t^{(1)})\\ + \frac{{\delta t}}{2}(1 - \frac{{{s_7}}}{2}){\partial _{1x}}(uF_{tx}^{(1)} - vF_{ty}^{(1)}) + \frac{{\delta t}}{2}(1 - \frac{{{s_8}}}{2}){\partial _{1y}}(uF_{ty}^{(1)} + vF_{tx}^{(1)}) = 0, \end{split} \end{equation}$ $

(44)

$\begin{equation} \begin{split} {\partial _{t2}}(\rho v) + \frac{1}{6}(1 - \frac{{{s_1}}}{2}){\partial _{1y}}{e^{(1)}} + (1 - \frac{{{s_7}}}{2})({\partial _{1x}}p_{xy}^{(1)} - \frac{1}{2}{\partial _{1y}}p_{xx}^{(1)}) \\ + \frac{{\delta t}}{2}(1 - \frac{{{s_7}}}{2}){\partial _{1x}}(uF_{ty}^{(1)} + vF_{tx}^{(1)})\\ + \frac{{\delta t}}{2}(1 - \frac{{{s_1}}}{2}){\partial _{1y}}(\mathbf{u} \cdot \mathbf{F}_t^{(1)}) - \frac{{\delta t}}{2}(1 - \frac{{{s_7}}}{2}){\partial _{1y}}(uF_{tx}^{(1)} - vF_{ty}^{(1)}) = 0. \end{split} \end{equation}$ $

(45)

To close the hydrodynamic equations at the second order of $ \varepsilon $, the terms of $ e^{(1)} $, $ p_{xx}^{(1)} $ and $ p_{xy}^{(1)} $ in Eqs. (44) and (45) should be estimated. Under the low Mach number assumption, these terms can be evaluated as

$\begin{equation} {e^{(1)}} = - \frac{1}{{{\widetilde{s}_1}}}[2\rho ({\partial _{1x}}u + {\partial _{1y}}v) + 3{s_1}\mathbf{u} \cdot \mathbf{F}_t^{(1)}] + O(M{a^3}), \end{equation}$ $

(46)

$\begin{equation} p_{xx}^{(1)} = - \frac{1}{{{\widetilde{s}_7}}}[\frac{2}{3}\rho ({\partial _{1x}}u - {\partial _{1y}}v) + {s_7}(uF_{tx}^{(1)} - vF_{ty}^{(1)})] + O(M{a^3}), \end{equation}$ $

(47)

$\begin{equation} p_{xy}^{(1)} = - \frac{1}{{{\widetilde{s}_8}}}[\frac{1}{3}\rho ({\partial _{1x}}v + {\partial _{1y}}u) + \frac{{{s_8}}}{2}(uF_{ty}^{(1)} + vF_{tx}^{(1)})] + O(M{a^3}). \end{equation}$ $

(48)

With Eqs. (46), (47) and (48), we can obtain the hydrodynamic equations at $ t_1 $ and $ t_2 $ scales,

Continuity equations

$\begin{equation} {\partial _{{t_1}}}\rho + {\partial _{1x}}(\rho {u}) + {\partial _{1y}}(\rho {v}) = 0, \end{equation}$ $

(49)

$\begin{equation} {\partial _{{t_2}}}\rho = 0. \end{equation}$ $

(50)

Momentum equations

$\begin{equation} {\partial _{{t_1}}}\rho {u} + {\partial _{1x}}(c_s^2\rho + \rho u^2) + {\partial _{1y}}(\rho {u}{v}) = F_{tx}^{(1)}, \end{equation}$ $

(51)

$\begin{equation} {\partial _{{t_1}}}\rho {v} + {\partial _{1x}}(\rho {u}{v}) + {\partial _{1y}}(c_s^2\rho + \rho v^2) = F_{ty}^{(1)}, \end{equation}$ $

(52)

$\begin{equation} {\partial _{t2}}(\rho u) = {\partial _{1x}}(\rho \nu ({\partial _{1x}}u - {\partial _{1y}}v) + \rho \zeta ({\partial _{1x}}u + {\partial _{1y}}v)) + {\partial _{1y}}(\rho \nu ({\partial _{1x}}v + {\partial _{1y}}u)), \end{equation}$ $

(53)

$\begin{equation} {\partial _{t2}}(\rho v) = {\partial _{1x}}(\rho \nu ({\partial _{1x}}v + {\partial _{1y}}u)) + {\partial _{1y}}(\rho \nu ({\partial _{1y}}v - {\partial _{1x}}u) + \rho \zeta ({\partial _{1x}}u + {\partial _{1y}}v)), \end{equation}$ $

(54)

where $ \nu = \rho \eta $ and $ \zeta = \rho \xi $ are the kinematic and bulk viscosities, respectively. In the present MRT-LB model, we enforce $ \nu = \frac{1}{3}(\frac{1}{s_7}-\frac{1}{2})\delta t $ and $ \zeta = \frac{1}{3}(\frac{1}{s_1}-\frac{1}{2})\delta t $.

Combining the above equations on $ t_0 $ and $ t_1 $ scale, the following hydrodynamic equations can be obtained.

$\begin{equation} {\partial _t}\rho + \nabla \cdot \left( {\rho {\mathbf{u}}} \right) = 0, \end{equation}$ $

(55)

$\begin{equation} \frac{{\partial \rho \mathbf{u}}}{{\partial t}} + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) = - \nabla c_s^2\rho + \nabla \cdot [\rho \nu (\nabla \mathbf{u} + \nabla {\mathbf{u}^T}) + \rho (\zeta - \nu )(\nabla \cdot \mathbf{u})\mathbf{I}] + {\mathbf{F}_t}. \end{equation}$ $

(56)

References

[1]	Clark W, Shih A, Hardin C, et al. (2003) Fixed abrasive diamond wire machining—part I: process monitoring and wire tension force. Int J Mach Tool Manu 43: 523–532. doi: 10.1016/S0890-6955(02)00215-8
[2]	Yu X, Wang P, Li X, et al. (2012) Thin Czochralski silicon solar cells based on diamond wire sawing technology. Sol Energ Mater Sol C 98: 337–342. doi: 10.1016/j.solmat.2011.11.028
[3]	Möller H (2006) Wafering of silicon crystals. Phys Stat Sol 203: 659–669. doi: 10.1002/pssa.200564508
[4]	Sopori B, Devayajanam S, Shet S, et al. (2013) Characterizing damage on Si wafer surfaces cut by slurry and diamond wire sawing. Proceedings of 39th IEEE PVSC 2013, Tampa, Florida, pp. 0945–0950.
[5]	Sopori B, Basnyat P, Devayajanam S, et al. (2015) Analyses of diamond wire sawn wafers: Effect of various cutting parameters. Proceedings of 42nd IEEE PVSC 2015, New Orleans, Louisiana.
[6]	De Meyer C, Heim B, Riddle Y (2012) Diamond wire wafering: A tutorial, 22^nd Silicon Workshop, July 22–25, 2012, Vail, Colorado, USA.
[7]	Bidiville A, Wasmer K, Kraft R, et al. (2009) Diamond wire-sawn silicon wafers – from the lab to the cell production. Presented at the 24th European Photovoltaic Solar Energy Conference and Exhibition, 21–25 September 2009, Hamburg, Germany, pp. 1400–1405.
[8]	Teomete E (2008) Mechanics of wire saw machining process: experimental analyses and modeling, Ph. D. thesis, Iowa State University.
[9]	Sopori B, Devayajanam S, Basnyat P (2015) Using minority carrier lifetime measurement to determine saw damage characteristics on Si wafer surfaces. Proceedings of 42nd IEEE PVSC 2015, New Orleans, Louisiana.
[10]	Wurzner S, Buchwald R, Mӧller H (2015) Surface damage and mechanical strength of silicon wafers. Phys Status Solidi C 12: 1119–1122. doi: 10.1002/pssc.201400227
[11]	Arif M, Rahman M, San W (2012) A state-of-the-art review of ductile cutting of silicon wafers for semiconductor and microelectronics industries. Int J Adv Manuf Technol 63:481–504. doi: 10.1007/s00170-012-3937-2
[12]	Blake P, Scattergood R (1990) Ductile- regime machining of germanium and silicon. J Am Ceram Soc 73: 949–957. doi: 10.1111/j.1151-2916.1990.tb05142.x
[13]	Ravindra D, Virkar S, Patten J (2011) Ductile mode micro laser assisted machining of silicon carbide in properties and applications of silicon carbide, Edited by Rosario Gerhardt, ISBN 978-953-307-201-2, 546 pages, Publisher: InTech, Chapters published.
[14]	Bidiville A, Wasmer K, Meer M, et al. (2015) Wire-sawing processes: parametrical study and modeling. Sol Energ Mater Sol C 132: 392–402. doi: 10.1016/j.solmat.2014.09.019
[15]	Donmich V, Gogotsi Y (2002) Phase transformation in silicon under contact loading. Rev Adv Mater Sci 3: 1–36. doi: 10.1016/S1468-6996(01)00150-4
[16]	Sopori B, Devayajanam S, Basnyat P (2016) A method for determining average damage depth of sawn crystalline silicon wafers. Rev Sc Inst 87: 45104. doi: 10.1063/1.4944792
[17]	Sopori B, Devayajanam S, Basnyat P, et al. (2015) Surface damage introduced by diamond wire sawing of Si wafers: Measuring in-depth and the lateral distributions for different cutting parameters. Mater Res Soc Symp Proc 1770: 61–66. doi: 10.1557/opl.2015.830
[18]	Choi C, Lee J, Cho S, et al. (1998) Evaluation of mechanical damage by high resolution X-ray diffraction and minority carrier lifetime. J Appl Phys 84: 168–173. doi: 10.1063/1.368073
[19]	Watanabe N, Kondo Y, Idle D, et al. (2010) Characterization of polycrystalline silicon wafers for solar cells sliced with novel fixed-abrasive wire. Prog Photovolt Res Appl 18: 485–490. doi: 10.1002/pip.923
[20]	Sopori B (1980) Rapid nondestructive technique for monitoring polishing damage in semiconductor wafers. Rev Sc Inst 50: 1513–1515.
[21]	Park H, Kwon S, Lee J, et al. (2009) Improvement on surface texturing of single crystalline silicon for solar cells by saw-damage etching using an acidic solution. Sol Energ Mater Sol C 93: 1773–1778.

This article has been cited by:

1.	Lei Xu, Zhengzheng Yan, Rongliang Chen, A discrete unified gas kinetic scheme on unstructured grids for viscid compressible flows and its parallel algorithm, 2023, 8, 2473-6988, 8829, 10.3934/math.2023443
2.	Chengjie Zhan, Zhenhua Chai, Baochang Shi, A lattice Boltzmann model for the coupled cross-diffusion-fluid system, 2021, 400, 00963003, 126105, 10.1016/j.amc.2021.126105
3.	Ben Gao, Qinglian Yin, Construction of invariant solutions and conservation laws to the $(2+1)$-dimensional integrable coupling of the KdV equation, 2020, 2020, 1687-2770, 10.1186/s13661-020-01466-6

Reader Comments

Your name:*

Email:*
© 2016 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)