Impact of particle shape on electron transport and lifetime in zinc oxide nanorod-based dye-sensitized solar cells

Roger Chang; Kemakorn Ithisuphalap; Ilona Kretzschmar; Roger Chang; Kemakorn Ithisuphalap; Ilona Kretzschmar

doi:10.3934/matersci.2016.1.51

AIMS Materials Science

2016, Volume 3, Issue 1: 51-65. doi: 10.3934/matersci.2016.1.51

Previous Article Next Article

Research article Topical Sections

Impact of particle shape on electron transport and lifetime in zinc oxide nanorod-based dye-sensitized solar cells

Department of Chemical Engineering, City College of the City University of New York, New York, NY 10031, USA

Received: 04 November 2015 Accepted: 03 January 2016 Published: 11 January 2016

Owing to its high electron mobility, zinc oxide represents a promising alternative to titanium dioxide as the working electrode material in dye-sensitized solar cells (DSCs). When zinc oxide is grown into 1-D nanowire arrays and incorporated into the working electrode of DSCs, enhanced electron dynamics and even a decoupling of electron transport (τ_d) and electron lifetime (τ_n) have been observed. In this work, DSCs with working electrodes composed of solution-grown, unarrayed ZnO nanorods are investigated. In order to determine whether such devices give rise to similar decoupling, intensity modulated photocurrent and photovoltage spectroscopies are used to measure τ_dand τ_n, while varying the illumination intensity. In addition, ZnO nanorod-based DSCs are compared with ZnO nanoparticle-based DSCs and nanomaterial shape is shown to affect electron dynamics. Nanorod-based DSCs exhibit shorter electron transport times, longer electron lifetimes, and a higher τ_n/τ_d ratio than nanoparticle-based DSCs.

Keywords:

Citation: Roger Chang, Kemakorn Ithisuphalap, Ilona Kretzschmar. Impact of particle shape on electron transport and lifetime in zinc oxide nanorod-based dye-sensitized solar cells[J]. AIMS Materials Science, 2016, 3(1): 51-65. doi: 10.3934/matersci.2016.1.51

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]	O'Regan B, Gratzel M (1991) A low-cost, high-efficiency solar cell based on dye-sensitized colloidal TiO₂ films. Nature 353: 734–740.
[2]	Mathew S, Yella A, Gao P, et al. (2014) Dye-sensitized solar cells with 13% efficiency achieved through the molecular engineering of porphyrin sensitizers. Nat Chem 6: 242–247. doi: 10.1038/nchem.1861
[3]	Kojima A, Teshima K, Shirai Y, et al. (2009) Organometal halide perovskites as visible-light sensitizers for photovoltaic cells. J Am Chem Soc 131: 6050–6051. doi: 10.1021/ja809598r
[4]	Lee MM, Teuscher J, Miyasaka T, et al. (2012) Efficient hybrid solar cells based on meso-superstructured organometal halide perovskites. Science 338: 643–647. doi: 10.1126/science.1228604
[5]	Burschka J, Pellet N, Moon S-J, et al. (2013) Sequential deposition as a route to high-performance perovskite-sensitized solar cells. Nature 499: 316–320. doi: 10.1038/nature12340
[6]	Yang WS, Noh JH, Jeon NJ, et al. (2015) High-performance photovoltaic perovskite layers fabricated through intramolecular exchange. Science 348: 1234–1237. doi: 10.1126/science.aaa9272
[7]	Xu F, Sun L (2011) Solution-derived ZnO nanostructures for photoanodes of dye-sensitized solar cells. Energy Environ Sci 4: 818–841. doi: 10.1039/C0EE00448K
[8]	Zhang Q, Cao G (2011) Nanostructured photoelectrodes for dye-sensitized solar cells. Nano Today 6: 91–109. doi: 10.1016/j.nantod.2010.12.007
[9]	Strehlow WH, Cook EL (1973) Compilation of energy band gaps in elemental and binary compound semiconductors and insulators. J Phys Chem Ref Data 2: 163–200. doi: 10.1063/1.3253115
[10]	Scanlon DO, Dunnill CW, Buckeridge J, et al. (2013) Band alignment of rutile and anatase TiO₂. Nat Mater 12: 798–801. doi: 10.1038/nmat3697
[11]	Chou TP, Zhang Q, Fryxell GE, et al. (2007) Hierarchically structured ZnO film for dye-sensitized solar cells with enhanced energy conversion efficiency. Adv Mater 19: 2588–2592. doi: 10.1002/adma.200602927
[12]	Dittrich T, Lebedev EA, Weidmann J (1998) Electron drift mobility in porous TiO₂ (anatase). Phys Status Solidi A 165: R5–R6.
[13]	Redmond G, Fitzmaurice D, Graetzel M (1994) Visible light sensitization by cis-bis(thiocyanato)bis(2,2'-bipyridyl-4,4'-dicarboxylato)ruthenium(II) of a transparent nanocrystalline ZnO film prepared by sol-gel techniques. Chem Mater 6: 686–691.
[14]	Hoyer P, Weller H (1995) Potential-dependent electron injection in nanoporous colloidal ZnO films. J Phys Chem 99: 14096–14100. doi: 10.1021/j100038a048
[15]	Rensmo H, Keis K, Lindström H, et al. (1997) High light-to-energy conversion efficiencies for solar cells based on nanostructured ZnO electrodes. J Phys Chem B 101: 2598–2601. doi: 10.1021/jp962918b
[16]	van de Lagemaat J, Frank AJ (2001) Nonthermalized electron transport in dye-sensitized nanocrystalline TiO₂ films: Transient photocurrent and random-walk modeling studies. J Phys Chem B 105: 11194–11205. doi: 10.1021/jp0118468
[17]	Frank A, Kopidakis N, Lagemaat Jvd (2004) Electrons in nanostructured TiO₂ solar cells: transport, recombination and photovoltaic properties. Coord Chem Rev 248: 1165–1179. doi: 10.1016/j.ccr.2004.03.015
[18]	Willis RL, Olson C, O'Regan B, et al. (2002) Electron dynamics in nanocrystalline ZnO and TiO₂ films probed by potential step chronoamperometry and transient absorption spectroscopy. J Phys Chem B 106: 7605–7613. doi: 10.1021/jp020231n
[19]	Quintana M, Edvinsson T, Hagfeldt A, et al. (2007) Comparison of dye-sensitized ZnO and TiO₂ solar cells: Studies of charge transport and carrier lifetime. J Phys Chem C 111: 1035–1041. doi: 10.1021/jp065948f
[20]	Kopidakis N, Benkstein KD, van de Lagemaat J, et al. (2003) Transport-limited recombination of photocarriers in dye-sensitized nanocrystalline TiO₂ solar cells. J Phys Chem B 107: 11307–11315.
[21]	Martinson ABF, McGarrah JE, Parpia MOK, et al. (2006) Dynamics of charge transport and recombination in ZnO nanorod array dye-sensitized solar cells. Phys Chem Chem Phys 8: 4655–4659. doi: 10.1039/b610566a
[22]	Baxter JB, Walker AM, Ommering Kv, et al. (2006) Synthesis and characterization of ZnO nanowires and their integration into dye-sensitized solar cells. Nanotechnol 17: S304. doi: 10.1088/0957-4484/17/11/S13
[23]	Galoppini E, Rochford J, Chen H, et al. (2006), Fast electron transport in metal organic vapor deposition grown dye-sensitized ZnO nanorod solar cells. J Phys Chem B 110: 16159–16161.
[24]	Law M, Greene LE, Johnson JC, et al. (2005) Nanowire dye-sensitized solar cells. Nat Mater 4: 455–459. doi: 10.1038/nmat1387
[25]	Noack V, Weller H, Eychmüller A (2002) Electron transport in particulate ZnO electrodes: A simple approach. J Phys Chem B 106: 8514–8523. doi: 10.1021/jp0200270
[26]	Qiu J, Li X, Zhuge F, et al. (2010) Solution-derived 40 µm vertically aligned ZnO nanowire arrays as photoelectrodes in dye-sensitized solar cells. Nanotechnol 21: 195602. doi: 10.1088/0957-4484/21/19/195602
[27]	Ku C-H, Wu J-J (2007) Chemical bath deposition of ZnO nanowire–nanoparticle composite electrodes for use in dye-sensitized solar cells. Nanotechnol 18: 505706. doi: 10.1088/0957-4484/18/50/505706
[28]	Xu F, Dai M, Lu Y, et al. (2010) Hierarchical ZnO nanowire−nanosheet architectures for high power conversion efficiency in dye-sensitized solar cells. J Phys Chem C 114: 2776–2782. doi: 10.1021/jp910363w
[29]	Xu C, Wu J, Desai UV, et al. (2011) Multilayer assembly of nanowire arrays for dye-sensitized solar cells. J Am Chem Soc 133: 8122–8125. doi: 10.1021/ja202135n
[30]	Seow ZLS, Wong ASW, Thavasi V, et al. (2009) Controlled synthesis and application of ZnO nanoparticles, nanorods and nanospheres in dye-sensitized solar cells. Nanotechnol 20: 045604. doi: 10.1088/0957-4484/20/4/045604
[31]	Cakir AC, Erten-Ela S (2012) Comparison between synthesis techniques to obtain ZnO nanorods and its effect on dye sensitized solar cells. Adv Powder Technol 23: 655–660. doi: 10.1016/j.apt.2011.08.003
[32]	Hosni M, Kusumawati Y, Farhat S, et al. (2014), Effects of oxide nanoparticle size and shape on electronic structure, charge transport, and recombination in dye-sensitized solar cell photoelectrodes. J Phys Chem C 118: 16791–16798.
[33]	Cao HL, Qian XF, Gong Q, et al. (2006) Shape- and size-controlled synthesis of nanometre ZnO from a simple solution route at room temperature. Nanotechnol 17: 3632. doi: 10.1088/0957-4484/17/15/002
[34]	Bai S, Hu J, Li D, et al. (2011) Quantum-sized ZnO nanoparticles: Synthesis, characterization and sensing properties for NO2. J Mater Chem 21: 12288–12294. doi: 10.1039/c1jm11302j
[35]	Hu X, Masuda Y, Ohji T, et al. (2008) Micropatterning of ZnO nanoarrays by forced hydrolysis of anhydrous zinc acetate. Langmuir 24: 7614–7617. doi: 10.1021/la8006348
[36]	Tan ST, Sun XW, Zhang XH, et al. (2006) Cluster coarsening in zinc oxide thin films by postgrowth annealing. J Appl Phys 100: 033502. doi: 10.1063/1.2218468
[37]	Hill JJ, Haller K, Gelfand B, et al. (2010) Eliminating capillary coalescence of nanowire arrays with applied electric fields. ACS Appl Mater Interfaces 2: 1992–1998. doi: 10.1021/am100290z
[38]	Sutton LE, HJM Bowen (1958) Tables of Interatomic Distances and Configuration in Molecules and Ions, London: The Chemical Society.
[39]	Spiro M, Creeth AM (1990) Tracer diffusion coefficients of I^-, I₃^-, Fe²⁺ and Fe³⁺ at low temperatures. J Chem Soc Faraday Trans 86: 3573–3576. doi: 10.1039/ft9908603573
[40]	Lin C, Tsai FY, Lee FH, et al. (2009) Enhanced performance of dye-sensitized solar cells by an Al₂O₃ charge-recombination barrier formed by low-temperature atomic layer deposition. J Mater Chem 19: 2999–3003. doi: 10.1039/b819337a
[41]	Mozer AJ, Wagner P, Officer DL, et al. (2008) The origin of open circuit voltage of porphyrin-sensitised TiO₂ solar cells. Chem Commun 4741–4743.
[42]	Ko SB, Cho AN, Kim MJ, et al., (2012) Alkyloxy substituted organic dyes for high voltage dye-sensitized solar cell: Effect of alkyloxy chain length on open-circuit voltage. Dyes Pigments 94: 88–98. doi: 10.1016/j.dyepig.2011.10.014
[43]	Kalyanasundaram K (2009) Dye-sensitized solar cells. Boca Raton, Fla: CRC.

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)