#  BAST: Binaural Audio Spectrogram Transformer for Binaural Sound
Localization

Sheng Kuang<EMAIL_ADDRESS>Kiki van der Heijden
<EMAIL_ADDRESS>Siamak Mehrkanoon
siamak.mehrkanoon<EMAIL_ADDRESS>Information and
Computing Sciences, Utrecht University, Utrecht, The Netherlands  Department
of Data Science and Knowledge Engineering, Maastricht University, The
Netherlands.  Donders Institute for Brain Cognition and Behavior, Radboud
University, Nijmegen, The Netherlands

###### Abstract

Accurate sound localization in a reverberation environment is essential for
human auditory perception. Recently, Convolutional Neural Networks (CNNs) have
been utilized to model the binaural human auditory pathway. However, CNN shows
barriers in capturing the global acoustic features. To address this issue, we
propose a novel end-to-end Binaural Audio Spectrogram Transformer (BAST) model
to predict the sound azimuth in both anechoic and reverberation environments.
Two modes of implementation, i.e. BAST-SP and BAST-NSP corresponding to BAST
model with shared and non-shared parameters respectively, are explored. Our
model with subtraction interaural integration and hybrid loss achieves an
angular distance of 1.29 degrees and a Mean Square Error of 1e-3 at all
azimuths, significantly surpassing CNN based model. The exploratory analysis
of the BAST’s performance on the left-right hemifields and anechoic and
reverberation environments shows its generalization ability as well as the
feasibility of binaural Transformers in sound localization. Furthermore, the
analysis of the attention maps is provided to give additional insights on the
interpretation of the localization process in a natural reverberant
environment.

###### keywords:

Transformer , Sound localization , Binaural integration

††journal: arXiv

##  1 Introduction

Sound source localization is a fundamental ability in everyday life.
Accurately and precisely localizing incoming auditory streams is required for
auditory perception and social communication. In the past decades, the
biological basis and neural mechanism of sound localization have been
extensively explored [1, 2, 3, 4]. Normal hearing listeners extract spatial
acoustic cues by mainly relying on the interaural level differences (ILD) and
interaural time differences (ITD) of the auditory input. These cues are
encoded through the human auditory subcortical pathway, in which the olivary
complex and associated physiological structures integrate and convey the
binaural signals from cochleas to the auditory cortex [2, 4]. However, sound
localization is frequently affected by the noise and reverberations in the
complex real-word environment, such as at a cocktail party. These noisy
signals may undercut or distort the spatial cues [5]. Although many studies
have focused largely on speech extraction from background noise, it is still
not clear how the spatial position of acoustic signals in complex listening
environments is extracted by the human brain.

Recently, Deep Learning (DL) [6] has been proposed to model auditory
processing and has achieved great success. These approaches enable
optimization of auditory models for real-life auditory environment [7, 8, 9,
10]. In the early attempts, DL methods were combined with conventional feature
engineering to deal with the noise and reverberation problems [7, 11]. For
instance, in [7], binaural spectral and spatial features were separately
extracted, providing complementary information for a two-layer Deep Neural
Network (DNN). Similarly, in [8, 12, 13], deep neural networks were used to
de-noise and de-reverberate complex sound stimuli. In a CNN-based azimuth
estimation approach, researchers utilized a Cascade of Asymmetric Resonators
with Fast-Acting Compression to analyze sound signals and used onsite-
generated correlograms to eliminate the echo interference [14]. However, most
of these approaches highly depend on feature selection. To reduce this
constraint, end-to-end Deep Residual Network (DRN) was recommended [15, 16].
Instead of selecting features from acoustic signal, raw spectrograms of sound
was utilized in Deep Residual Network for azimuth prediction [16]. DRN was
shown to be robust even in the presence of unknown noise interference at the
low signal-to-noise ratio. Subsequently, [9] proposed a pure attention-based
Audio Spectrogram Transformer (AST) and achieved the state-of-the-art results
for audio classification on multiple datasets. Although these DL-based methods
have yielded promising results, however due to a lack of similar architectures
to the human binaural auditory pathway, they may not resemble the neural
processing of sound localization.

To encode the neural mechanisms underlying sound localization, the performance
of deep learning methods is commonly compared to the human behavioral
characteristics [17, 18, 19]. For instance, [19] systematically explored the
performance of binaural sound clips localization of a CNN in a real-life
listening environment, however, its architecture does not resemble the
structure of human auditory pathway. This issue has been addressed by
utilizing a hierarchical neurobiological-inspired CNN (NI-CNN) to model the
binaural characteristics of human spatial hearing [17]. This unique
hierarchical design, models the binaural signal integration process and is
shown to have brain-like latent feature representations. However, NI-CNN [17]
is not an end-to-end model as it leverages a cochlear method to generate
auditory nerve representations as model input. Furthermore, considering the
wide range of frequencies of sound input, the convolution operations in NI-CNN
mainly extract local-scale features and therefore may exhibit limitations for
extracting global features in the acoustic time-frequency spectrogram.

In this study, we build on the success and barriers of previously proposed
deep neural networks at localizing sound sources to further develop an end-to-
end transformer based model, which captures the global acoustic features from
auditory spectrograms. We aimed at (i) investigating the performance of a pure
Transformer-based hierarchical binaural neural network for addressing real-
life sound localization; (ii) exploring the effect of various loss functions
and binaural integration methods on the localization acuity at different
azimuths; (iii) visualizing the attention flow of the proposed model to
demonstrate the localization process.

##  2 Related Works

Binaural auditory models utilize head-related binaural sound waves to
accomplish specific acoustic tasks. Gaussian Mixture Model (GMM), commonly
used in machine learning-based studies, calculates the probability
distribution of the source location in reverberant environments [20].
Subsequently, model-based methods have been encouraged to extract ILD and ITD
cues for DNN training [7]. However, the performance of these hybrid techniques
remains unstable since the feature extraction routine varies across different
datasets. NI-CNN, as a data-driven DL method, is able to learn latent features
for azimuth prediction from human auditory nerve representations [17].
However, it still requires a computational neuroscience-based cochlear model
to generate these representations. In this paper, we design an end-to-end
model to automatically learn features from spectrogram without using model-
based methods.

Transformer was originally proposed in natural language processing to handle
long-range dependencies [21, 22]. Recently, Transformer was successfully
applied in the computer vision field by casting images into patch embedding
sequences [23, 24]. Many hybrid models combined the Transformer with a CNN or
Recurrent Neural Network (RNN) in audio processing and some studies even
directly embedded attention blocks into CNN or RNN to capture global features
in a parameter-efficient way [25, 26, 27].

Subsequently, the authors in [9] introduced AST which uses Transformer model
and variable-length monaural spectrograms to perform sound classification
tasks. As a convolution–free, pure attention-based model, AST uses the
overlapped-patch embedding generation policy to convert the intra-patch local
features to inter-patch attention weights. AST has achieved the state-of-the-
art results on multiple datasets [28, 29] for audio classification tasks.

Figure 1: Architecture of the proposed Binaural Audio Spectrogram Transformer
(BAST). (a) The architecture of the proposed model. (Here there are
$N_{H}N_{T}$ number of patches). (b) The architecture of a single Transformer
encoder. (c) Three examined interaural integration methods: concatenation,
addition and subtraction.

##  3 Method

###  3.1 Model architecture

The proposed Binaural Audio Spectrogram Transformer (BAST), is illustrated in
Fig. 1. Similar to NI-CNN [17], a dual-input hierarchical architecture is
utilized to simulate the human subcortical auditory pathway. As opposed to NI-
CNN which uses convolution layers, here three Transformer Encoders (i.e.,
left, right and center), hereafter called TE-L, TE-R and TE-C are utilized to
construct a pure attention-based model. In particular, the pre-processed left
and right sound waves are converted to left and right spectrograms denoted by
$x^{L}\in\mathbb{R}^{H\times T}$ and $x^{R}\in\mathbb{R}^{H\times T}$. Here,
$H$ indicates the frequency band and $T$ indicates the number of Tukey windows
( with shape parameter: 0.25).

In what follows, the TE-L path to process the input data is explained. The
other path, i.e TE-R, follows the same process. At the beginning of patch
embedding layer, the left spectrogram $x^{L}\in\mathbb{R}^{H\times T}$ is
first decomposed into an overlapped-patch sequence
$x_{patch}^{L}\in\mathbb{R}^{P^{2}\times(N_{H}N_{T})}$, where $P$ is the patch
size, $N_{H}$ and $N_{T}$ are the number of patches in height and width
respectively obtained as follows,

$N_{H}=\left\lceil\frac{H-P+S}{S}\right\rceil,N_{T}=\left\lceil\frac{T-P+S}{S}\right\rceil.$
(1)

In case $H-P$ and $T-P$ are not divisible by the stride $S$ between patches,
the spectrogram is zero-padded on the top and right respectively. A trainable
linear projection is added to flatten each patch to a $D$ dimensional latent
representation, hereafter called patch embeddings [23]. Since our model
outputs the sound location coordinates, the classification token in the
Transformer encoder is removed. A fixed absolute position embedding [23] is
added to the patch embeddings to capture the position information of the
spectrogram in the Transformer. Here, the learnable position embedding is not
used as it did not significantly change model performance compared to absolute
position embedding[9]. The output of the left position embedding layer
$z_{in}^{L}\in\mathbb{R}^{D\times(N_{H}N_{T})}$ is then fed to the Transformer
encoder TE-L.

We use the identical Transformer encoder design in [23, 9], consisting of $K$
stacked Multi-head Self-Attention (MSA) and Multi-Layer Perceptron (MLP)
blocks. The BAST model performance is compared when using shared and non-
shared parameters across the left and right Transformer encoders. Hereafter,
BAST-SP refers to BAST model whose left and right Transformer encoders share
parameters whereas in BAST-NSP the parameters of the left and right
Transformer encoders are not shared. The output of left and right Transformer
encoder, i.e., $z_{out}^{L}$ and $z_{out}^{R}$, represents neural signals
underlying the initial auditory processing stage along the left and right
auditory pathway respectively. Subsequently, these binaural feature maps are
integrated to simulate the function of the human olivary nucleus. Similar to
NI-CNN, here three integration methods, i.e. addition, subtraction and
concatenation, are investigated. Specifically, addition is the summation of
feature maps of both sides; subtraction represents left feature map subtracted
from the right feature map; concatenation is implemented by concatenating
$z_{out}^{L}$ and $z_{out}^{R}$ along the the first dimension to produce
$z_{in}^{C}\in\mathbb{R}^{2D\times(N_{H}N_{T})}$. TE-C receives the integrated
feature map $z_{in}^{C}$ and output sequence $z_{out}^{C}$. Next, an average
operation of the patch dimension and a linear transformer layer is applied to
finally produce the sound location coordinates $(x,y)$. The last linear layer
does not have any activation function, therefore the estimated coordinates can
be any point on the 2D plane.

###  3.2 Loss Function

Three loss functions, i.e., Angular Distance (AD) loss [30], Mean Square Error
(MSE) loss, as well as hybrid loss with a convex combination of AD and MSE,
are explored in training the proposed model. Let $C_{i}=(x_{i},y_{i})$ and
$\hat{C}_{i}=(\hat{x}_{i},\hat{y}_{i})$ denote the ground truth and
predication coordinates for the i-$th$ sample. MSE loss measures the squared
difference of Euclidean distance between the prediction and the ground truth
as follows:

$\textrm{MSE}=\frac{1}{N}\sum_{i}^{N}\|C_{i}-\hat{C}_{i}\|_{2}^{2},$ (2)

where $(\hat{x_{i}},\hat{y_{i}})$ is the predicted coordinates,
$(x_{i},y_{i})$ is the true sound location and $N$ is the batch size. Note
that MSE loss is able to penalize the large Euclidean distance error but is
insensitive to the angular distance, which means that the azimuth may differ
with the same MSE. In contrast to MSE, AD loss merely measures the angular
distance while ignoring the Euclidean distance:

$\textrm{AD}=\frac{1}{\pi
N}\sum^{N}_{i}\arccos{(\frac{C_{i}\hat{C}_{i}^{T}}{\|C_{i}\|_{2}\|\hat{C}_{i}\|_{2}})},$
(3)

where $C_{i},\hat{C}_{i}\neq 0$.

Table 1: The performance of BAST-NSP and BAST-SP compared to NI-CNN and NI-CNN∗ when different loss and binaural integration methods are used. The best performed model in AD and MSE are in bold. The ↓ indicates the lower the value of the metric, the better the model performance. Model | Loss | Angular Distance(AD) ↓ | Mean Squared Error (MSE) ↓
---|---|---|---
Concatenation | Addition | Subtraction | Concatenation | Addition | Subtraction
NI-CNN[17] | MSE | 4.80° | 4.80° | 5.30° | 0.011 | 0.013 | 0.014
AD | 3.70° | 3.90° | 5.20° | — | — | —
NI-CNN[17]∗ | MSE | 8.92° | 3.51° | 3.67° | 0.077 | 0.032 | 0.038
AD | 7.85° | 1.97° | 1.85° | — | — | —
Hybrid | 8.35° | 3.53° | 3.19° | 0.074 | 0.033 | 0.031
BAST-NSP | MSE | 2.78° | 2.48° | 2.42° | 0.003 | 0.002 | 0.002
AD | 2.39° | 1.30° | 1.63° | — | — | —
Hybrid | 2.76° | 1.83° | 1.29° | 0.004 | 0.002 | 0.001
BAST-SP | MSE | 2.02° | 4.97° | 1.94° | 0.002 | 0.018 | 0.002
AD | 2.66° | 13.87° | 1.43° | — | — | —
Hybrid | 1.98° | 5.72° | 2.03° | 0.003 | 0.026 | 0.002
Table 2: The number of layers in each Transform encoder as well as the total number of trainable parameters of the proposed models. Tuple ($\cdot$, $\cdot$, $\cdot$) indicates the number of layers in the left, right and center Transformer encoder respectively. Model |  | Interaural
---
Integration
| Transformer
---
Layers
| Trainable
---
Parameters
BAST-NSP | Concatenation | (3, 3, 3) | ~76M
| Addition/
---
Subtraction
(3, 3, 3) | ~57M
BAST-SP | Concatenation | (3, 3, 3) | ~57M
| Addition/
---
Subtraction
(3, 3, 3) | ~38M

##  4 Experiments

###  4.1 Dataset

We use the binaural audio data in [17], which consists of a training dataset
and an independent testing dataset. In the training dataset, 4600 real-life
sound waves (duration:500 ms, sampling rate: 16000) are placed in 36 azimuth
positions, respectively, with $10\degree$ azimuth resolution, $0\degree$
elevation, and 1-meter distance from the center point. In addition, sound
waves are spatialized with two acoustic environments, i.e. an anechoic
environment (AE) without reverberation and a 10m $\times$ 14m lecture hall
with reverberation (RV). In particular, here the training and test sets
contain data from both AE and RV environments. In total, the training dataset
has 331200 binaural learning samples. Similarly, the independent testing
dataset contains 400 new sound waves processed with the same method as
described above, producing 28800 testing samples.

###  4.2 Model Evaluation

Models were evaluated by means of MSE and AD errors defined in Eq. (2) and
(3). The lower AD and MSE errors the better localization performance is. Note
that the MSE metric has no meaning when BAST is trained by AD loss because
this loss does not optimize the Euclidean distance between the ground truth
and the prediction, and that BAST has no constraints on the numerical range of
the predicted coordinates.

###  4.3 Training Settings

As mentioned in 3.1, each sound wave are transformed to binaural spectrogram
(size: 2$\times$129$\times$61, frequency range: 0-8000Hz, window length:
128ms, overlap: 64ms) before training. In order to have balanced training
samples, we randomly select 75% binaural spectrograms in each azimuth position
and listening environments of training dataset. The remaining data is used as
validation set. As stated before in the Dataset section, a separate test
dataset is available for this study. This setting results in $248400$ training
samples and $82800$ validation samples. Here, Adam optimizer [31] is used to
train the model for 50 epochs with a batch size of 48 and a fixed learning
rate of 1e-4. In patch embedding layers, the stride of patches is set to 6,
yielding 180 patches per spectrogram. Each Transformer encoder has three
layers, with 1024 hidden dimensions (2048 dimensions when using concatenation
as integration method in the last Transformer encoder TE-C), 16 attention
heads in MSA blocks, 1024 dimensions in MLP blocks, 0.2 dropout rate in patch
embeddings and MLP blocks. Our model implementation 111code available at
https://github.com/ShengKuangCN/BAST is based on Python 3.8 and Pytorch 1.9.0,
and are trained from scratch on 2 $\times$ NVIDIA GeForce GTX 1080Ti GPUs with
11GB of RAM. The empirically found and used number of layers in each
Transformer encoder as well as the total number of trainable parameters are
presented in Table 2.

Figure 2: The angular distance (AD) error of the proposed BAST in each azimuth
with different loss functions and interaural integration methods.

##  5 Results

###  5.1 Overall Performance

The performance of the proposed BAST-NSP and BAST-SP models are compared with
those of NI-CNN model. In particular, for NI-CNN model two modes of
implementations corresponding to cochleogram and spectrogram as model inputs
have been considered and are denoted by NI-CNN and NI-CNN∗ respectively. The
obtained results of the compared models with different combinations of
binaural integration methods and loss functions are tabulated in Table 1.
BAST-NSP has achieved the best AD error of 1.29°and the best MSE of 0.001 when
using the subtraction binaural integration and the hybrid loss function.
Compared to NI-CNN, BAST-NSP reduces AD error 65.4% from 3.70°to 1.29°and MSE
90.9% from 0.011 to 0.001. In addition, BAST-NSP outperforms NI-CNN∗ (NI-CNN∗:
AD=1.85°, MSE=0.031), although both models received the same input. BAST-SP
has achieved AD=1.43°, MSE=0.002, therefore surpassing the performance of NI-
CNN and NI-CNN∗ models while it is still inferior to BAST-NSP. Overall, in
this study, BAST-NSP outperforms the other tested models in performing
binaural sound localization.

We further analyze the influence of different binaural integration methods on
the BAST-NSP and BAST-SP performance. Here, the performances are compared in
terms of AD error. Specifically, in both cases when the BAST-NSP is trained by
AD loss and hybrid loss, binaural integration through addition and subtraction
improved the model performance compared to concatenation (AD loss: Add.=1.30°,
Sub.=1.63°, Concat.=2.39°; Hybrid loss: Add.=1.83°, Sub. =1.29°,
Concat.=2.76°, see Table 1). In case of MSE loss, the performance across the
three integration methods of BAST-NSP is similar. In BAST-SP, addition
integration causes a huge AD error increment over BAST-NSP (BAST-SP: 4.97°,
BAST-NSP: 1.30°), indicating that the left-right identical feature addition
brings a great challenge to the model to predict the azimuth.

The effect of three types of loss functions on BAST-NSP and BAST-SP
performance are not the same. In BAST-NSP, AD loss achieves lower AD when
using concatenation or addition, but the hybrid loss yields the lowest AD in
terms of subtraction. In BAST-SP, one can observe an interaction of loss
function and binaural integration methods, i.e. the best loss function depends
on the applied binaural integration methods.

###  5.2 Performance at different azimuths

To better understand the localization performance of the models in different
azimuth, the test AD error of each azimuth is shown in Fig. 2. The test AD
error in BAST-NSP is much smaller when the sound source is located closer to
the interaural midline. However, this pattern is not observed in BAST-SP.

Figure 3: The AD error of the proposed BAST-NSP and BAST-SP in the left and
right hemifield. The boxplot indicates quartiles of the metric distribution
with respect to azimuths. The asterisk between two boxes indicates the
statistical significance (p$<$0.05, paired t-test with FDR correction) between
the left and right hemifield. Figure 4: The MSE of the proposed BAST-NSP and
BAST-SP in the left and right hemifield. The boxplot indicates quartiles of
the metric distribution with respect to azimuths. The asterisk between two
boxes indicates the statistical significance (p$<$0.05, paired t-test with FDR
correction) between the left and right hemifield.  Figure 5: An example of the
attention matrices in the proposed model (i.e., BAST-NSP, hybrid loss and
subtraction). The corresponding sound clip was randomly selected in the
category of human speech with reverberation. For each layer, we present the
patch-to-patch attention matrix (size: 180$\times$180) calculated by the
rollout method in [32]. Note that we initialize the attention matrix at the
first layer of TE-C by summing the attention matrices at the last layer of
TE-L and TE-R. Figure 6: Attention rollout corresponding to the spectrogram
shown in Fig. 5. (a): Left and right spectrogram. (b): The left and right
attention rollout obtained from the 3rd layer of TE-L and TE-R Transformers.
(c): The left and right Attention rollout obtained from the 3rd layer of TE-C.

###  5.3 Performance in left and right hemifield

To explore the symmetry of the model predictions, we further compare the
evaluation metrics between the left-right hemifields. One can observe
comparable model performance in left and right hemifield in Fig. 3 and 4. This
result is confirmed by paired t-test (False Discovery Rate (FDR) corrected for
multiple comparisons). More specifically, Fig. 3, shows an insignificant
difference of AD error between the left and right hemifield in most conditions
(corrected p$>$0.05). However, a minor but significant difference (corrected
p$<$0.05) is observed in BAST-NSP trained with MSE loss and addition
integration, and in BAST-SP trained with AD loss and subtraction. More
precisely, the difference in AD error between the left and right hemifield was
not significant in most conditions (corrected p$>$0.05, Fig. 3), thus
supporting the consistent symmetry of model predictions.

###  5.4 Performance in different environments

We conduct two additional experiments to illustrate the generalization of the
proposed model by training in one listening environment and testing in both
environments separately, i.e., AE and RV. This analysis is conducted on the
best performing model, i.e., BAST-NSP with hybrid loss and subtraction
integration method. As shown in Table 3, the model that is trained using the
data of both AE and RV environments, achieves the best test results compared
to other models which are trained only using the data of one of the
environment.

Table 3: The performance of the proposed BAST-NSP model in different listening environments. AE and RV indicate the anechoic and reverberation environments respectively.  | Training
---
Environment
| Testing
---
Environment
AD | MSE
AE | AE | 1.14° | 0.001
RV | 8.66° | 0.027
RV | AE | 16.70° | 0.078
RV | 1.65° | 0.002
AE+RV | AE | 1.10° | 0.001
RV | 1.48° | 0.001

###  5.5 Attention Analysis

To interpret the localization process, we utilize Attention Rollout [33] to
visualize the attention maps of the proposed model (BAST-NSP with subtraction
method and hybrid loss). Rollout calculates the attention matrix by
recursively multiplying the attention matrices along the forward propagation
path. [32] enhanced this method by adding an additional identical matrix
before multiplication to simulate the effect of residual connection of MSA.
Due to the interaural integration layer, in BAST-NSP, we initialize the
attention matrix of TE-C by summing the attention weights from both sides
regardless of the integration method.

Fig. 5 shows the patch-to-patch attention matrices (size: 180$\times$180) of a
randomly selected spectrogram. In the first layer of each Transformer, most of
the patches are self-focused and pay attention to some scattered patches.
However, in the last layer, all patches yield nearly consistent attention
weights to some specific patches. The attention rollout heat map with respect
to the left and right spectrogram are depicted in Fig. 6. Although parameter-
sharing setting is not used in BAST-NSP, one can still observe that the model
focuses most of its attention on similar regions on both sides, see Fig. 6(b).
The final attention map, Fig. 6(c), shows that the model further processes the
attention after the integration layer and boosts the attention weights in
bottom left regions.

##  6 Conclusion

In this paper, a novel Binaural Audio Spectrogram Transformer (BAST) for sound
source localization is proposed. The obtained results show that this pure
attention-based model leads to significant azimuth acuity improvement compared
to NI-CNN model. In particular, subtraction interaural integration and hybrid
loss is the best training combination for BAST. Additionally, we found that
the performance and statistical significance in left-right hemifields vary
with different combinations of training settings. In conclusion, this work
contributes to a convolution-free model of real-life sound localization. The
data and implementation of our BAST model are available at
https://github.com/ShengKuangCN/BAST.

## References

  * [1] D. W. Batteau, The role of the pinna in human localization, Proceedings of the Royal Society of London. Series B. Biological Sciences 168 (1011) (1967) 158–180.
  * [2] J. O. Pickles, Auditory pathways: anatomy and physiology, Handbook of clinical neurology 129 (2015) 3–25.
  * [3] K. van der Heijden, J. P. Rauschecker, B. de Gelder, E. Formisano, Cortical mechanisms of spatial hearing, Nature Reviews Neuroscience 20 (10) (2019) 609–623.
  * [4] B. Grothe, M. Pecka, D. McAlpine, Mechanisms of sound localization in mammals, Physiological reviews 90 (3) (2010) 983–1012.
  * [5] J. Blauert, S. Hearing, The psychophysics of human sound localization, in: Spatial Hearing, MIT Press, 1997.
  * [6] Y. LeCun, Y. Bengio, G. Hinton, Deep learning, nature 521 (7553) (2015) 436–444.
  * [7] X. Zhang, D. Wang, Deep learning based binaural speech separation in reverberant environments, IEEE/ACM transactions on audio, speech, and language processing 25 (5) (2017) 1075–1084.
  * [8] S. Y. Lee, J. Chang, S. Lee, Deep learning-based method for multiple sound source localization with high resolution and accuracy, Mechanical Systems and Signal Processing 161 (2021) 107959.
  * [9] Y. Gong, Y.-A. Chung, J. Glass, Ast: Audio spectrogram transformer, arXiv preprint arXiv:2104.01778 (2021).
  * [10] T.-D. Truong, C. N. Duong, H. A. Pham, B. Raj, N. Le, K. Luu, et al., The right to talk: An audio-visual transformer approach, in: Proceedings of the IEEE/CVF International Conference on Computer Vision, 2021, pp. 1105–1114.
  * [11] L. Perotin, R. Serizel, E. Vincent, A. Guérin, Crnn-based joint azimuth and elevation localization with the ambisonics intensity vector, in: 2018 16th International Workshop on Acoustic Signal Enhancement (IWAENC), IEEE, 2018, pp. 241–245.
  * [12] T. Yoshioka, S. Karita, T. Nakatani, Far-field speech recognition using cnn-dnn-hmm with convolution in time, in: 2015 IEEE international conference on acoustics, speech and signal processing (ICASSP), IEEE, 2015, pp. 4360–4364.
  * [13] S. Park, Y. Jeong, H. S. Kim, Multiresolution cnn for reverberant speech recognition, in: 2017 20th Conference of the Oriental Chapter of the International Coordinating Committee on Speech Databases and Speech I/O Systems and Assessment (O-COCOSDA), IEEE, 2017, pp. 1–4.
  * [14] Y. Xu, S. Afshar, R. K. Singh, R. Wang, A. van Schaik, T. J. Hamilton, A binaural sound localization system using deep convolutional neural networks, in: 2019 IEEE International Symposium on Circuits and Systems (ISCAS), IEEE, 2019, pp. 1–5.
  * [15] K. He, X. Zhang, S. Ren, J. Sun, Deep residual learning for image recognition, in: Proceedings of the IEEE conference on computer vision and pattern recognition, 2016, pp. 770–778.
  * [16] N. Yalta, K. Nakadai, T. Ogata, Sound source localization using deep learning models, Journal of Robotics and Mechatronics 29 (1) (2017) 37–48.
  * [17] K. van der Heijden, S. Mehrkanoon, Goal-driven, neurobiological-inspired convolutional neural network models of human spatial hearing, Neurocomputing 470 (2022) 432–442.
  * [18] K. van der Heijden, S. Mehrkanoon, Modelling human sound localization with deep neural networks., in: ESANN, 2020, pp. 521–526.
  * [19] A. Francl, J. H. McDermott, Deep neural network models of sound localization reveal how perception is adapted to real-world environments, Nature Human Behaviour 6 (1) (2022) 111–133.
  * [20] N. Ma, J. A. Gonzalez, G. J. Brown, Robust binaural localization of a target sound source by combining spectral source models and deep neural networks, IEEE/ACM Transactions on Audio, Speech, and Language Processing 26 (11) (2018) 2122–2131.
  * [21] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, I. Polosukhin, Attention is all you need, Advances in neural information processing systems 30 (2017).
  * [22] J. Devlin, M.-W. Chang, K. Lee, K. Toutanova, Bert: Pre-training of deep bidirectional transformers for language understanding, arXiv preprint arXiv:1810.04805 (2018).
  * [23] A. Dosovitskiy, L. Beyer, A. Kolesnikov, D. Weissenborn, X. Zhai, T. Unterthiner, M. Dehghani, M. Minderer, G. Heigold, S. Gelly, et al., An image is worth 16x16 words: Transformers for image recognition at scale, arXiv preprint arXiv:2010.11929 (2020).
  * [24] K. Han, Y. Wang, H. Chen, X. Chen, J. Guo, Z. Liu, Y. Tang, A. Xiao, C. Xu, Y. Xu, et al., A survey on vision transformer, IEEE Transactions on Pattern Analysis and Machine Intelligence (2022).
  * [25] Z. Zhang, S. Xu, S. Zhang, T. Qiao, S. Cao, Attention based convolutional recurrent neural network for environmental sound classification, Neurocomputing 453 (2021) 896–903.
  * [26] Y.-B. Lin, Y.-C. F. Wang, Audiovisual transformer with instance attention for audio-visual event localization, in: Proceedings of the Asian Conference on Computer Vision, 2020.
  * [27] Q. Kong, Y. Xu, W. Wang, M. D. Plumbley, Sound event detection of weakly labelled data with cnn-transformer and automatic threshold optimization, IEEE/ACM Transactions on Audio, Speech, and Language Processing 28 (2020) 2450–2460.
  * [28] K. J. Piczak, Esc: Dataset for environmental sound classification, in: Proceedings of the 23rd ACM international conference on Multimedia, 2015, pp. 1015–1018.
  * [29] P. Warden, Speech commands: A dataset for limited-vocabulary speech recognition, arXiv preprint arXiv:1804.03209 (2018).
  * [30] X. Xiao, S. Zhao, X. Zhong, D. L. Jones, E. S. Chng, H. Li, A learning-based approach to direction of arrival estimation in noisy and reverberant environments, in: 2015 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), IEEE, 2015, pp. 2814–2818.
  * [31] D. P. Kingma, J. Ba, Adam: A method for stochastic optimization, arXiv preprint arXiv:1412.6980 (2014).
  * [32] H. Chefer, S. Gur, L. Wolf, Transformer interpretability beyond attention visualization, in: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2021, pp. 782–791.
  * [33] S. Abnar, W. Zuidema, Quantifying attention flow in transformers, arXiv preprint arXiv:2005.00928 (2020).


# A diffuse-interface lattice Boltzmann method for fluid-particle interaction
problems

Jiao Liu  Changsheng Huang  Zhenhua Chai<EMAIL_ADDRESS>Baochang Shi
School of Mathematics and Statistics, Huazhong University of Science and
Technology, Wuhan 430074, China  School of Mathematics and Physics, China
University of Geosciences, Wuhan 430074, China  Hubei Key Laboratory of
Engineering Modeling and Scientific Computing, Huazhong University of Science
and Technology, Wuhan 430074, China

###### Abstract

In this paper, a diffuse-interface lattice Boltzmann method (DI-LBM) is
developed for fluid-particle interaction problems. In this method, the sharp
interface between the fluid and solid is replaced by a thin but nonzero
thickness transition region named diffuse interface, where the physical
variables varies continuously. In order to describe the diffuse interface, we
introduce a smooth function, which is similar to the order parameter in phase-
field model or the volume fraction of solid phase in the partially saturated
lattice Boltzmann method (PS-LBM). In addition, to depict the fluid-particle
interaction more accurately, a modified force term is also proposed and
included in the evolution equation of the DI-LBM. Some classical problems are
used to test the DI-LBM, and the results are in good agreement with some
available theoretical and numerical works. Finally, it is also found that the
DI-LBM is more efficient and accurate than the PS-LBM with the superposition
model.

###### keywords:

diffuse-interface lattice Boltzmann method , fluid-particle interaction ,
Chapman-Enskog analysis

##  1 Introduction

Particulate flows are ubiquitous in both nature and engineering, such as, the
sediment deposition, fluidized beds, and so on [1, 2]. With the rapid
development of computer technology, the numerical simulation has been becoming
an important and efficient tool in the study of the particulate flows. The
lattice Boltzmann method (LBM), as a kinetic-based numerical approach, has
gained a great success in the simulation of complex hydrodynamic problems [3,
4, 5, 6]. Compared to the traditional methods for the Navier-Stokes equations,
the LBM has some distinct advantages, including the clear physical background,
easy implementation of boundary conditions, simplicity in programming and high
computational efficiency [3]. Due to these advantages of the LBM, it has also
been applied to investigate particulate flows [7, 8, 9]. Generally speaking,
in the framework of LBM, there are three basic methods in the treatment of
fluid-particle interface, namely, the bounce-back method, the immersed
boundary (IB) method and the partially saturated (PS) method.

Ladd [8, 9] first proposed a shell method to simulate fluid-particle flows
based on the half-way bounce-back scheme. In the bounce-back scheme, the
inside and outside of the particle are filled with fluid, and the same bounce-
back operation is performed for both the internal and external fluids of the
particle. However, this method requires the particle boundary to be in the
middle of the grid, which may cause the simulated particle boundary to be
different from the actual physical boundary. To overcome this problem, several
interpolation-based curved boundary schemes have been developed and applied to
describe the particle-fluid interaction [11, 12, 13]. In this kind of method,
there is no fluid inside the particle, and it has a second-order accuracy for
arbitrary curved boundaries. However, the hydrodynamic force at the boundary
is usually not smooth, and the fictitious oscillation is more serious when the
moving particle is considered [14].

The IB method was first developed by Peskin [15] to study blood flows in the
heart. Later, Feng [16] proposed an immersed-boundary lattice Boltzmann method
(IB-LBM) to simulate fluid-particle interaction problems. In this method, the
Euler mesh is used for fluid domain, while the particle boundary is marked by
a set of Lagrangian points. The fluid-particle interaction is achieved by
adding an external force to the fluid with the Dirac function. However, the
IB-LBM usually dose not satisfy the no-slip boundary condition [19]. To
enforce the no-slip boundary at the fluid-particle interface, some implicit
IB-LBMs have also been developed [17, 18, 19].

The partially saturated (PS) method coupled with the LBM (hereafter the PS-
LBM) was proposed by Noble and Torczynski [20], and then it is also applied to
study some two-dimensional flows by Cook et al. [21]. In this method, the
fluid is filled with the whole domain. However, unlike above two methods, the
fluid-particle interaction in this method is realized through modifying the
collision term in the evolution equation of LBM where a parameter representing
the volume fraction of solid phase is introduced. According to the modified
collision term, there are two different models, i.e., the non-equilibrium
bounce back (BB) and the superposition (SP) models. Compared to the bounce-
back method [8, 9], this method can give a relatively smooth hydrodynamic
force.

Owing to the advantage in mass conservation and without special treatment on
the fresh fluid nodes [6], the PS-LBM has received increasing attention in
studying the particulate flows [23, 22]. However, it should be noted that
there is another problem in the PS-LBM, i.e., the accurate computation of the
solid-phase volume fraction is very complicated, especially for three-
dimensional problems. To eliminate this limitation of the PS-LBM, in this work
a smooth function is first introduced to represent the volume fraction of
solid phase, and then, based on the PS-LBM with SP model, we developed a
diffuse-interface LBM (DI-LBM) for fluid-particle interaction problems. The
rest of the paper is organized as follows. In Section 2, a diffuse-interface
lattice Boltzmann method is developed. In Section 3, some numerical examples
are carried out to test the present DI-LBM, and finally, some conclusions are
given in Section 4.

##  2 Numerical method

In this section, we first presented a brief introduction to the standard LBM
and PS-LBM, and then following the PS-LBM with SP model, a DI-LBM is developed
for fluid-particle interaction problems.

###  2.1 Standard lattice Boltzmann method

In the standard LBM for fluid flows, the evolution equation reads [3]

$f_{i}(\mathbf{x}+\textbf{c}_{i}\Delta t,t+\Delta
t)=f_{i}(\textbf{x},t)-{1\over\tau}[f_{i}(\mathbf{x},t)-{f_{i}}^{eq}(\textbf{x},t)],$
(1)

where $f_{i}(\textbf{x},t)$ is the distribution function at position x and
time $t$. $\\{\mathbf{c}_{i},i=0,\cdots,q-1\\}$ is the set of discrete
velocity, $q$ is the number of discrete velocity directions, $\Delta t$ is the
time step, $\tau$ is the dimensionless relaxation time.
$f_{i}^{eq}(\textbf{x},t)$ is the equilibrium distribution function, and is
defined by

$f_{i}^{eq}(\textbf{x},t)={\omega_{i}}\rho\left[{1+{{{\textbf{c}_{i}}\cdot\textbf{u}}\over{c_{s}^{2}}}+{({{\textbf{c}_{i}}\cdot\textbf{u}})^{2}\over{2c_{s}^{4}}}-{{{\textbf{u}^{2}}}\over{2c_{s}^{2}}}}\right],$
(2)

where $\omega_{i}$ is the weight coefficient, ${c_{s}}$ is the lattice sound
speed, and is related to the lattice speed $c$ ($c={{\Delta
x}\mathord{\left/{\vphantom{{\Delta x}{\Delta t}}}\right.\kern-1.2pt}{\Delta
t}}$, $\Delta x$ is the lattice spacing).

For simplicity but without loss of generality, in this work we only considered
the two-dimensional problems, and adopted the commonly used D2Q9 lattice model
[24], in which the discrete velocities are defined as

$\textbf{c}_{i}=\begin{cases}(0,0),&i=0,\\\
(\cos[(i-1)\pi/2],\sin[(i-1)\pi/2])c,&i=1,2,3,4,\\\
\sqrt{2}(\cos[(2i-1)\pi/4],\sin[(2i-1)\pi/4])c,&i=5,6,7,8.\end{cases}$ (3)

The corresponding weight coefficients are given by

$w_{i}=\begin{cases}4/9,&i=0,\\\ 1/9,&i=1,2,3,4,\\\
1/36,&i=5,6,7,8.\end{cases}$ (4)

In addition, the macroscopic variable $\rho$ and u can be computed by

$\rho=\sum_{i}f_{i},\;\rho\textbf{u}=\sum_{i}\textbf{c}_{i}f_{i}.$ (5)

Finally, through the Chapman-Enskog analysis [5], the Navier-Stokes equations
can be recovered from Eq. (1) with the following viscosity,

$\upsilon=c_{s}^{2}\left(\tau-{1\over 2}\right)\Delta t.$ (6)

###  2.2 The partially saturated lattice Boltzmann method

In the PS-LBM, the fluid-particle interaction is realized by introducing an
additional collision term $\Omega_{i}^{s}$, and consequently, the evolution
equation can be written as

$f_{i}(\mathbf{x}+\textbf{c}_{i}\Delta t,t+\Delta
t)=f_{i}(\textbf{x},t)+(1-B)\left\\{-{1\over\tau}[f_{i}(\mathbf{x},t)-{f_{i}}^{eq}(\textbf{x},t)]\right\\}+B\Omega_{i}^{s},$
(7)

where $B$ is a weighting function of the solid-phase volume fraction
$\varepsilon_{s}$, and usually it can be simply given by

$B(\varepsilon_{s},\tau)=\varepsilon_{s}.$ (8)

It can be find that $\varepsilon_{s}=0$ ($B=0$) denotes the pure fluid, and in
this case, Eq. (7) would degenerate into the standard LBM (1); while
$\varepsilon_{s}=1$ ($B=1$) represents the pure solid. For the SP model in the
PS-LBM, the additional collision term $\Omega_{i}^{s}$ in Eq. (7) is designed
as

$\Omega_{i}^{s}=f_{i}^{eq}(\rho,\textbf{u}_{s})-f_{i}(\textbf{x},t)+\big{(}1-{1\over\tau}\big{)}[f_{i}(\textbf{x},t)-f_{i}^{eq}(\rho,\textbf{u})],$
(9)

where $\textbf{u}_{s}$ is the velocity of the particle.

The total hydrodynamic force and torque imposed on solid particle are
calculated by summing the additional collision operator over the nodes covered
by the solid,

$\textbf{F}_{f}=-{\Delta x^{2}\over\Delta
t}\sum_{n}B_{n}\sum_{i}\Omega_{i}^{s}\textbf{c}_{i},$ (10)
$\textbf{T}_{f}=-{\Delta x^{2}\over\Delta
t}\sum_{n}(\textbf{x}_{n}-\textbf{x}_{s})\times\left(B_{n}\sum_{i}\Omega_{i}^{s}\textbf{c}_{i}\right),$
(11)

where $\textbf{x}_{n}$ is the coordinate of the node $n$, and $\textbf{x}_{s}$
is the position of the particle. In addition, we would also like to point out
that Eq. (7) can also be written into another form [22],

$f_{i}(\mathbf{x}+\textbf{c}_{i}\Delta t,t+\Delta
t)=f_{i}(\textbf{x},t)-{1\over\tau}[f_{i}(\mathbf{x},t)-{f_{i}}^{eq}(\textbf{x},t)]+B[f_{i}^{eq}(\rho,\textbf{u}_{s})-f_{i}^{eq}(\rho,\textbf{u})].$
(12)

It is clear that compared to the evolution equation (1) in the standard LBM,
only an additional term
$B[f_{i}^{eq}(\rho,\textbf{u}_{s})-f_{i}^{eq}(\rho,\textbf{u})]$ is included
in Eq. (12).

###  2.3 A diffuse-interface lattice Boltzmann method

From above discussion, one can find that the PS-LBM remains two advantages of
the standard LBM, i.e., the locality of the collision operator and simple
linear streaming operator, and thus it is easy to solve the problems involving
the moving boundaries [25]. However, the accurate calculation of local solid-
phase volume fraction $\varepsilon_{s}$ in the PS-LBM is complicated, and the
computational cost is also expensive. To overcome these problems, here we
introduce a parameter $\phi$, similar to the order parameter in the diffuse-
interface method for multiphase flows [26] or the solid-phase volume fraction
$\varepsilon_{s}$ in the PS-LBM, and develop a new DI-LBM for particulate
flows. Similar to the PS-LBM, in the DI-LBM, the sharp fluid-particle
interface is replaced by a thin but nonzero transition region where the
physical variables change continuously. In addition, based on the PS-LBM with
the SP model [20] or Eq. (12), the evolution of the DI-LBM is written as

$f_{i}(\mathbf{x}+\textbf{c}_{i}\Delta t,t+\Delta
t)=f_{i}(\textbf{x},t)+\left\\{-{1\over\tau}[f_{i}(\mathbf{x},t)-{f_{i}}^{eq}(\textbf{x},t)]\right\\}+\phi
F_{i},$ (13)

where the parameter $\phi$ is a hyperbolic tangent function, and is defined as

$\phi={1+\tanh(2l/\varepsilon)\over 2},$ (14)

where $l$ is the distance to the boundary, $\varepsilon$ is the thickness of
diffuse interface. $F_{i}$ is the discrete force term used to replace the last
term $[f_{i}^{eq}(\rho,\textbf{u}_{s})-f_{i}^{eq}(\rho,\textbf{u})]$ in Eq.
(12), and is given by

$F_{i}=\big{(}1-{1\over
2\tau}\big{)}w_{i}\rho\left({\textbf{c}_{i}\cdot(\textbf{u}_{s}-\textbf{u}^{*})\over
c_{s}^{2}}+{(\textbf{u}_{s}\textbf{u}_{s}-\textbf{u}^{*}\textbf{u}^{*}):(\textbf{c}_{i}\textbf{c}_{i}-c_{s}^{2}\textbf{I})\over
2c_{s}^{4}}\right),$ (15)

where $\textbf{u}^{*}$ is the velocity without considering the fluid-particle
interaction, u is the corrected velocity, and they are defined by

$\textbf{u}^{*}={\sum_{i}\textbf{c}_{i}f_{i}\over\rho},$ (16)
$\textbf{u}=\textbf{u}^{*}+{1\over 2}\phi(\textbf{u}_{s}-\textbf{u}^{*}).$
(17)

The hydrodynamic force $\textbf{F}_{f}$ and torque $\textbf{T}_{f}$ can be
calculated through the first-order moment of the discrete force term,

$\textbf{F}_{f}=-{\Delta x^{2}\over\Delta
t}\sum_{n}\phi_{n}\sum_{i}F_{i}\textbf{c}_{i},$ (18) $\textbf{T}_{f}=-{\Delta
x^{2}\over\Delta
t}\sum_{n}(\textbf{x}_{n}-\textbf{x}_{s})\times\left(\phi_{n}\sum_{i}F_{i}\textbf{c}_{i}\right).$
(19)

The calculation procedure of the DI-LBM is as follows:

  1. 1.

Initialize the values of variables.

  2. 2.

Compute the distribution functions in entire computational domain using Eq.
(13).

  3. 3.

Correct the velocity field using Eq. (17).

  4. 4.

Calculate the hydrodynamic force and torque according to Eqs. (18) and (19).

  5. 5.

Update the velocity of the boundary and the parameter $\phi$.

  6. 6.

Repeat steps 2-5 until convergence is reached.

It should be noted that in the DI-LBM, a smooth function $\phi$ is used to
replace the volume fraction of solid phase in the PS-LBM, and the difficulty
in the computation of the solid-phase volume fraction can be avoided.
Therefore, compared to the original PS-LBM with the SP model, the present DI-
LBM would be more efficient. In addition, we proposed a modified discrete
force term [see Eq. (15)], which also makes the present LB method more
accurate than the original PS-LBM in the study of the particulate flows (see
the numerical results in the following section). Through the Chapman-Enskog
analysis (see the details in the Appendix), we can obtain the macroscopic
equations of the present DI-LBM, which are similar to those of the IB-LBM
[19].

##  3 Numerical results and discussion

In this section, we conducted some numerical tests to validate the present DI-
LBM; including the fluid flows passing a stationary circular cylinder, a
particle settling along channel centerline, the off-centerline sedimentation
of a particle, the sedimentation of two particles and a neutrally buoyant
particle moving in the Poiseuille flow. In the following simulations, the non-
equilibrium extrapolation scheme [27] is applied to treat the physical
boundary conditions of above problems.

###  3.1 The fluid flows passing a stationary circular cylinder

We first considered the problem of fluid flows passing a stationary circular
cylinder, and the schematic of the problem is shown in Fig. 1. It is known
that the flow behavior of the problem is mainly governed by the Reynolds
number ($Re$), which is a dimensionless number defined by

$Re={UD\over\upsilon},$ (20)

where $U$ is the free stream velocity, $D$ is the diameter of the cylinder,
$\upsilon$ is the kinematic viscosity. In addition, for a specified case where
the density $\rho$, the diameter $D$ and the velocity $U$ are given, we can
calculate the drag coefficient,

$C_{D}={F_{fx}\over 0.5\rho DU^{2}},$ (21)

where $F_{fx}$ is the $x$-component of the hydrodynamic force.

Figure 1: The schematic of fluid flows passing a stationary circular cylinder.

In our simulations, the computational domain is $35D\times 20D$, the circular
cylinder with the diameter $D=1.0$ is placed at $(10D,\,10D)$, the free stream
velocity is set to be $U=0.1$. We performed some simulations, and plotted the
streamlines in Fig. 2 where $Re=20$ and $40$. From this figure, it can be
observed that when the flow field reaches steady state, there is a pair of
symmetric recirculating eddies formed behind the cylinder, and the length of
the recirculating region increases with the increase of $Re$. We noted that
these results are qualitatively consistent with the previous work [28]. In
addition, to give a quantitative comparison between the present results and
some available data [28, 29], we also calculated the drag coefficient and
dimensionless recirculation length $L_{w}=2L/D$ ($L$ is the recirculation
length), and presented them in Table 1. It can be seen from this table that
the results of the present DI-LBM are in good agreement with those reported in
the previous studies [28, 29].

Figure 2: The streamlines of the fluid flows around the stationary circular cylinder at $Re=20$ (a) and $40$ (b). Table 1: A comparison of the drag coefficient and recirculation length between the present work and some previous studies. $Re$ |  | References | $C_{D}$ | $2L/D$
---|---|---|---|---
$Re=20$ |  | Present | $2.172$ | $1.875$
|  | Niu et al. [28] | $2.144$ | $1.89$
|  | He et al. [29] | $2.152$ | $1.842$
$Re=40$ |  | Present | $1.64$ | $4.796$
|  | Niu et al. [28] | $1.589$ | $4.52$
|  | He et al. [29] | $1.499$ | $4.490$

###  3.2 A particle settling along channel centerline

We continued to study the problem of a particle settling along channel
centerline. The configuration of this problem is depicted in Fig. 3 where a
particle with the diameter $D$ is placed in the center of the channel, the
width of the channel is $W$. Initially, the particle is released from rest,
and then it would sink under the gravity force. It should be noted that unlike
the previous problem, besides the Navier-Stokes equations, we also need to
solve the following equations for particle motion,

$m_{s}{d\textbf{u}_{s}\over
dt}=\textbf{F}_{f}+\left({\rho_{f}\over\rho_{s}}-1\right)m_{s}\textbf{g},$
(22) $I_{s}{d\omega_{s}\over dt}=T_{f},$ (23) ${d\textbf{x}_{s}\over
dt}=\textbf{u}_{s},$ (24)

where $\textbf{u}_{s}$ is the velocity of the setting particle, $m_{s}$ and
$\omega_{s}$ are the mass and the angular velocity of the particle, g is the
gravity acceleration, $I_{s}$ is the rotational inertia of the particle.

Theoretically, for the problem with a small Reynolds number, one can obtain
the approximate solution of the steady drag force [30],

Figure 3: The configuration of a particle settling along channel centerline.

$\textbf{F}_{d}=4\pi K\mu\textbf{u}_{s}.$ (25)

The parameter $K$ is a correction factor, and is given by

$K=(\ln
W^{*}-0.9175+1.7244W^{*-2}-1.7302W^{*-4}+2.4056W^{*-6}-4.5913W^{*-8})^{-1},$
(26)

where $W^{*}=W/D$ is the blockage ratio.

When the problem reaches the equilibrium state, we can also derive the final
velocity of the settling particle [31],

$\textbf{u}_{s}={D^{2}\over 16K\mu}(\rho_{f}-\rho_{s})\textbf{g}.$ (27)

In our simulations, the diameter of the cylinder $D$ is set to be $0.24\,\rm
cm$, the computational domain is $W\times H=1.2\,\rm cm\times 6.0\,\rm cm$,
the density of the fluid is $\rho_{f}=1.0\,\rm g/\rm cm^{3}$, and the
viscosity of the fluid is $\mu=0.1\,\rm g/(\rm cm\cdot\rm s)$. At the initial
time $t=0$, the particle is released at the point $(0.6\,\rm cm,\,3.0\,\rm
cm)$. For the left and right channel walls, the no-slip boundary condition is
applied. At the inlet, the velocity is set to be zero, while at the outlet,
the fully developed boundary conditions is adopted.

Figure 4: The settling velocities of the particle at different values of the
particle density ($\rho_{s}=1.01,1.02$ and $1.05\rm g/\rm cm^{3}$).

We carried out some simulations with three different values of the particle
density ($\rho_{s}=1.01,1.02$ and $1.05\,\rm g/\rm cm^{3}$), and presented the
settling velocities of three different cases in Fig. 4 where the grid size is
$\Delta x=0.01$ and the relaxation time is $\tau=0.8$. As shown in this
figure, the velocity of the settling particle increases with the increase of
particle density, and the present results also agree well with the previous
numerical data [31]. In addition, it is also found that when the particle
density is small, the numerical results are very close to the theoretical
solution (27), while if the particle density becomes large (e.g.,
$\rho_{s}=1.05\,\rm g/\rm cm^{3}$), there would be a deviation between the
numerical and theoretical results. This is because with the increase of the
particle density, the settling velocity of the particle and the Reynolds
number would be increased, which may cause the theoretical prediction (27) to
be inaccurate.

###  3.3 The off-centerline sedimentation of a particle

Next, we considered the sedimentation of a particle off the centerline. The
schematic of the problem is shown in Fig. 5 where a particle with the diameter
$D=0.1\,\rm cm$ is located at the off-center line of the channel. The width
and height of the channel are $W=4D$ and $H=400/13D$, the initial position of
the particle is $(0.19W,\,0.75H)$, the density of the fluid is
$\rho_{f}=1.0\,\rm g/\rm cm^{3}$, and the viscosity of fluid is $0.01\,\rm
cm^{2}/s$. Similar to above discussion, here two different values of the
particle density, i.e., $\rho_{s}=1.01$ and $1.03\,\rm g/\rm cm^{3}$ are
considered.

Figure 5: The schematic of the off-centerline sedimentation of a particle
settling in a channel.

Figure 6:  The trajectory of the particle and time history of the particle
velocity [(a): trajectory; (b): $x$-component velocity; (c): $y$-component
velocity; (d): angular velocity].

We conducted some simulations with the lattice size $Nx\times Ny=104\times
800$ and the relaxation time $\tau=0.6$, and presented the trajectory of the
particle denoted by $(x_{s},y_{s})$ and time history of the particle velocity
in Fig. 6. Based on the results shown in Fig. 6(a), one can find that the
particle finally moves to the centerline of the channel $(x_{c}=0.2\,\rm cm)$,
which is due to the asymmetric force induced by the channel walls. From Fig.
6(b) and 6(d), we can also observe that with the increase of time, the
horizontal velocity $u_{s}$ and the angular velocity $\omega_{s}$ gradually
become zero. This is because when the time is large enough, the particle would
reach the steady state along the channel centerline. In addition, the results
in Fig. 6(c) indicate that with the increase of the particle density, the
vertical velocity $v_{s}$ increase owing to the large gravity. Finally, we
would like to point out that the present results agree well with the numerical
data reported in the previous work [32].

###  3.4 The sedimentation of two particles

To further test the capacity of the DI-LBM in the study of the particulate
flows, the problem of the sedimentation of two particles is also considered.
For this problem, the diameter of the particles is $D=0.2\,\rm cm$, the
channel size is $2\,\rm cm\times 8\,\rm cm$, the viscosity of fluid is
$0.01\,\rm g/(\rm cm\cdot\rm s)$, the densities of the fluid and particles are
$\rho_{f}=1.0\,\rm g/\rm cm^{3}$ and $\rho_{s}=1.01\,\rm g/\rm cm^{3}$. The
initial positions of two particles are $(0.999\,\rm cm,\,7.2\,\rm cm)$ and
$(1.0\,\rm cm,\,6.8\,\rm cm)$.

In the following simulations, to avoid the overlap between two particles or
between the particles and the walls, some short-range repulsive forces should
be included. In this work, we adopted the collision model proposed by Wan and
Turek [33], in which the repulsive forces can be expressed as

$F_{i,j}^{P-P}=\begin{cases}0,&d_{i,j}>r_{i}+r_{j}+\xi,\\\
{1\over\epsilon_{p}^{{}^{\prime}}}(X_{i}-X_{j})(r_{i}+r_{j}-d_{i,j}),&d_{i,j}\leq
r_{i}+r_{j},\\\ {1\over\epsilon_{p}}(X_{i}-X_{j})(r_{i}+r_{j}+\xi-
d_{i,j})^{2},&r_{i}+r_{j}\leq d_{i,j}\leq r_{i}+r_{j}+\xi,\end{cases}$ (28)

$F_{i}^{P-W}=\begin{cases}0,&d_{i}^{{}^{\prime}}>2r_{i}+\xi,\\\
{1\over\epsilon_{W}^{{}^{\prime}}}(X_{i}-X_{i}^{{}^{\prime}})(2r_{i}-d_{i}^{{}^{\prime}}),&d_{i}^{{}^{\prime}}\leq
2r_{i},\\\ {1\over\epsilon_{W}}(X_{i}-X_{i}^{{}^{\prime}})(2r_{i}+\xi-
d_{i}^{{}^{\prime}})^{2},&2r_{i}\leq d_{i}^{{}^{\prime}}\leq
2r_{i}+\xi,\end{cases}$ (29)

where $r_{i}$ is the radius of the $i$th particle, $d_{i,j}=|X_{i}-X_{j}|$,
$d_{i}^{{}^{\prime}}=|X_{i}-X_{i}^{{}^{\prime}}|$. $X_{i}$ is the $i$th
particle center, $X_{i}^{{}^{\prime}}$ is the coordinate vector of the center
of the nearest imaginary particle on the boundary wall. $\epsilon_{P}$ and
$\epsilon_{P}^{{}^{\prime}}$ are small positive stiffness parameters for
particle-particle collisions, and are fixed as $1.0\times 10^{-7}$.
$\epsilon_{W}$ and $\epsilon_{W}^{{}^{\prime}}$ are two small positive
stiffness parameters for particle-wall collisions, and are taken as
$\epsilon_{W}=\epsilon_{P}/2$ and
$\epsilon_{W}^{{}^{\prime}}=\epsilon_{P}^{{}^{\prime}}/2$. $\xi$ is the
threshold that is set to be one lattice unit.

Figure 7: The vertical positions of the two particles. Figure 8: The
horizontal positions of the two particles.

Figs. 7 and 8 show the vertical and horizontal positions of the two particles,
which are also compared with the results of Jafari et al. [34]. From these two
figures, one can find that as time goes on, the trailing particle first comes
close to the leading particle because of the low pressure wake created by the
leading particle, then the trailing particle with a larger velocity induces a
kissing contact with the leading particle, and finally two particles tumble
and separate from each other. These three distinct processes are so-called
drafting, kissing and tumbling. In addition, it is also observed that the
present results are close to those in the previous study [34].

###  3.5 A neutrally buoyant particle moving in the Poiseuille flow

The last problem we considered is the motion of a neutrally buoyant particle
in the Poiseuille flow. The configuration of this problem is depicted in Fig.
9 where computational domain is $L\times W=20\times 4$, the diameter of the
particle is $D=1$, and the initial position of the particle is
$(x_{0},y_{0})=(L/2,\,0.25W)$. In order to reduce the computational cost, the
moving computational domain is adopted, when the particle moves one lattice
unit in the horizontal direction, the computational domain also moves one
lattice in the same direction [31]. The pressure drop from inlet to outlet is
$\triangle p=P_{in}-P_{out}=0.00267$, the lattice spacing is $\Delta x=1/25$,
and the relaxation time is $\tau=0.75$.

Figure 9: Schematic of a particle moving in the Poiseuille flow. Figure 10: The trajectory of the moving particle in the Poiseuille flow. Table 2: A comparison of the equilibrium position of the particle between the present work and some previous studies.  Present |  | PS-LBM |  | Li et al. [32] |  | Tao et al. [14]
---|---|---|---|---|---|---
$0.2876$ |  | $0.50$ |  | $0.2874$ |  | $0.2849$

We carried out some simulations, and presented the trajectory of the particle
in Fig. 10. As shown in this figure, with the increase of time, the particle
would finally reach the equilibrium position between the wall and the
centerline of the channel, which is also in good agreement with the previous
work [14]. However, when the original PS-LBM with the SP model is adopted to
investigate this problem, the particle would reach the centerline of the
channel, which may be caused by the inappropriate collision term or the
discrete force term [20, 22]. In addition, we also measured the equilibrium
positions of the particle with the DI-LBM and PS-LBM with the SP model, and
conducted a comparison of the values obtained with different numerical methods
in Table 2. From this table, one can also observe that the present result
agrees well with the previous works [32, 14], while the original PS-LBM with
the SP model cannot give the accurate result.

##  4 Conclusions

In this work, we proposed a diffuse-interface lattice Boltzmann method for
fluid-particle interaction problems. The distinct feature of the method is
that the sharp boundary between the fluid and solid is represented by a
diffuse interface with a finite thickness, and simultaneously a smooth
function $\phi$ is introduced to describe the diffuse interface. This method
can not only overcome the disadvantage of the PS-LBM in the computation of
solid-phase volume fraction, but also improve the accuracy of the PS-LBM
through developing a modified force term. The present DI-LBM is validated by
several numerical experiments, and the results show that the DI-LBM is
effective and accurate in the study of the fluid-particle interaction
problems.

In addition, we would also like to point out that the present DI-LBM can be
considered as an effective tool for more complex fluid-particle interaction
problems, for instance, the moving particles in porous media [35].

## Acknowledgements

This work was financially supported by the National Natural Science Foundation
of China (Grants No. 12072127, No. 11702259 and No. 51836003).

## Appendix

We now perform a detailed Chapman-Enskog analysis to obtain the macroscopic
equations from the DI-LBM (13). In the Chapman-Enskog analysis, the
distribution function, the time and space derivatives, and the force term can
be first expanded as

$\displaystyle f_{i}=f_{i}^{(0)}+\epsilon
f_{i}^{(1)}+\epsilon^{2}f_{i}^{(2)}+\cdots,$ (A.1a)
$\displaystyle\partial_{t}=\epsilon\partial_{t_{1}}+\epsilon^{2}\partial_{t_{2}},\;\;\nabla=\epsilon\nabla_{1},\;\;F_{i}=\epsilon
F_{i}^{(1)},$ (A.1b)

where $\epsilon$ is a small parameter proportional to the Knudsen number.
Then, taking the Taylor expansion to Eq. (13) , we have

$\Delta tD_{i}f_{i}+{\Delta t^{2}\over
2}{D_{i}}^{2}f_{i}+\cdots=\left\\{-{1\over\tau}[f_{i}(\mathbf{x},t)-{f_{i}}^{eq}(\textbf{x},t)]\right\\}+\phi
F_{i}.$ (A.2)

If we substitute Eq. (A.1) into Eq. (A.2), one can obtain

$O(\epsilon^{0}):f_{i}^{(0)}=f_{i}^{eq}(\rho,\textbf{u}),$ (A.3)
$O(\epsilon^{1}):\quad D_{1i}f_{i}^{(0)}=-{1\over\tau\Delta
t}f_{i}^{(1)}+{\phi\over\Delta t}F_{i}^{(1)},$ (A.4)
$O(\epsilon^{2}):\quad\partial_{t_{2}}f_{i}^{(0)}+D_{1i}f_{i}^{(1)}+{\Delta
t\over 2}D_{1i}^{2}f_{i}^{(0)}=-{1\over\tau\Delta t}f_{i}^{(2)},$ (A.5)

where $D_{1i}=\partial_{t_{1}}+\mathbf{c}_{i}\cdot\nabla_{1}$. According to
Eq. (A.3), we can derive the following moments,

$\begin{split}&\sum_{i}f_{i}^{(0)}=\rho\\\
&\sum_{i}\mathbf{c}_{i}f_{i}^{(0)}=\rho\textbf{u},\\\
&\sum_{i}\mathbf{c}_{i}\mathbf{c}_{i}f_{i}^{(0)}=c_{s}^{2}\rho\mathbf{I}+\rho\textbf{u}\textbf{u},\\\
&\sum_{i}\mathbf{c}_{i}\mathbf{c}_{i}\mathbf{c}_{i}f_{i}^{(0)}=c_{s}^{2}\Delta\cdot\rho\textbf{u},\end{split}$
(A.6)

which can be used to obtain the zero and first-order moments of the non-
equilibrium distribution functions,

$\sum_{i}f_{i}^{(k)}=0\;(k\geq
1),\;\;\sum_{i}\mathbf{c}_{i}f_{i}^{(1)}=-{\phi\over
2}\rho(\textbf{u}_{s}-\textbf{u}^{*}),\sum_{i}\mathbf{c}_{i}f_{i}^{(2)}=0.$
(A.7)

Substituting Eq. (A.4) into Eq. (A.5) yields

$O(\epsilon^{2}):\quad\partial_{t_{2}}f_{i}^{(0)}+D_{1i}(1-{1\over
2\tau})f_{i}^{(1)}=-{1\over\tau\Delta t}f_{i}^{(2)}-D_{1i}{\phi\over
2}F_{i}^{(1)}.$ (A.8)

Summing Eqs. (A.4) and (A.8) over $i$, we have

$\partial_{t_{1}}\rho+\nabla_{1}\cdot(\rho\textbf{u})=0,$ (A.9)
$\partial_{t_{2}}\rho=0.$ (A.10)

Combining Eqs. (A.9) and (A.10), one can obtain the continuity equation,

$\partial_{t}\rho+\nabla\cdot(\rho\textbf{u})=0.$ (A.11)

If we multiply $\mathbf{c}_{i}$ on both sides of Eqs. (A.4) and (A.8), one can
derive

$\partial_{t_{1}}(\rho\textbf{u})+\nabla_{1}\cdot(c_{s}^{2}\rho\mathbf{I}+\rho\textbf{u}\textbf{u})=\bar{F}^{(1)},$
(A.12) $\partial_{t_{2}}(\rho\textbf{u})+\nabla_{1}\cdot\left(1-{1\over
2\tau}\right)\sum\mathbf{c}_{i}\mathbf{c}_{i}f_{i}^{(1)}=-\nabla_{1}\cdot{\phi\over
2}\sum\mathbf{c}_{i}\mathbf{c}_{i}F_{i}^{(1)},$ (A.13)

where $\bar{F}^{(1)}={\phi\over\Delta t}\rho(\textbf{u}_{s}-\textbf{u}^{*})$.

With the help of Eq. (A.4), we can give an evaluation to second-order moment
of $f_{i}^{(1)}$,

$\begin{split}-{1\over\tau\Delta
t}\sum\mathbf{c}_{i}\mathbf{c}_{i}f_{i}^{(1)}&=\sum\mathbf{c}_{i}\mathbf{c}_{i}D_{1i}f_{i}^{eq}-{\phi\over\Delta
t}\sum\mathbf{c}_{i}\mathbf{c}_{i}F_{i}^{(1)}\\\ &\approx
c_{s}^{2}\rho[\nabla_{1}\textbf{u}+(\nabla_{1}\textbf{u})^{T}].\\\
\end{split}$ (A.14)

Substituting Eq. (A.14) into Eq. (A.13), and combining the equations at
$t_{1}$ and $t_{2}$ scales, we have

$\begin{split}\partial_{t}(\rho\textbf{u})+\nabla\cdot(\rho\textbf{u}\textbf{u})&=-\nabla
P+\nabla\cdot\left\\{\big{(}\tau-{1\over 2}\big{)}\Delta
tc_{s}^{2}\rho\left[\nabla\textbf{u}+(\nabla\textbf{u})^{T}\right]\right\\}+\textbf{f},\end{split}$
(A.15)

where f is the force caused by the fluid-particle interaction.

From above discussion, one can find that the macroscopic equations (A.11) and
(A.15) can be obtained from the DI-LBM with the following viscosity,

$\mu=c_{s}^{2}\rho\left(\tau-{1\over 2}\right)\Delta t.$ (A.16)

## Acknowledgements

## References

  * [1] S. Li, J.S. Marshall, G. Liu, Q. Yao, Adhesive particulate flow: The discrete-element method and its application in energy and environmental engineering, Prog. Energy Combust. Sci. 37 (6) (2011) 633-668.
  * [2] M. Maxey, Simulation methods for particulate flows and concentrated suspensions, Annu. Rev. Fluid Mech. 49 (2017) 171-193.
  * [3] S.Y. Chen, G.D. Doolen, Lattice Boltzmann method for fluid flows, Annu. Rev. Fluid Mech. 30 (1998) 329-364.
  * [4] S. Succi, _The Lattice Boltzmann Equation for Fluid Dynamics and Beyond_ (Oxford University Press, Oxford, 2001).
  * [5] Z. Guo, C. Shu, _Lattice Boltzmann Method and Its Applications in Engineering_ (World Scientific Publishing Co. Pte. Ltd., Singapore, 2013).
  * [6] T. Krüger, H. Kusumaatmaja, G. Silva, O. Shardt, A. Kuzmin, E.M. Viggen, _The Lattice Boltzmann Method: Principles and Practice_ (Springer International Publishing, Switzerland, 2017).
  * [7] C. K. Aidun, J. R. Clausen, Lattice-Boltzmann method for complex flows, Annu. Rev. Fluid Mech. 42 (2010) 439-472.
  * [8] A.J.C. Ladd, Numerical simulations of particulate suspensions via a discretized Boltzmann equation, J. Fluid Mech. 271 (1994) 285-309, 311-339.
  * [9] A.J.C. Ladd, R. Verberg, Lattice-Boltzmann simulations of particle-fluid suspensions, J. Stat. Phys. 104 (2001) 1191-1251.
  * [10] C.K. Aidun, Y. Lu, Lattice Boltzmann simulation of solid particles suspended in fluid, J. Stat. Phys. 81 (1995) 49-61.
  * [11] D. Yu, R. Mei, L.S. Luo, W. Shyy, Viscous flow computations with the method of lattice Boltzmann equation, Progress in Aerospace Sciences 39 (2003) 329-367.
  * [12] M. Bouzidi, M. Firdaouss, P.J. Lallemand, Momentum transfer of a Boltzmann-lattice fluid with boundaries, J. Comput. Phys. 13 (2003) 3452-3459.
  * [13] W. Zhao, W.A. Yong, Single-node second-order boundary schemes for the lattice Boltzmann method, J. Comput. Phys. 329 (2017) 1-15.
  * [14] S. Tao, J.J. Hu, Z.L. Guo, An investigation on momentum exchange methods and refilling algorithms for lattice Boltzmann simulation of particulate flows, Comput. Fluids 133 (2016) 1-14.
  * [15] C.S. Peskin, Numerical analysis of blood flow in the heart, J. Comput. Phys. 25 (3) (1977) 220-252.
  * [16] Z.G. Feng, E.E. Michaelides, The immersed boundary-lattice Boltzmann method for solving fluid-particles interaction problems, J. Comput. Phys. 195 (2) (2004) 602-628.
  * [17] J. Wu, C. Shu, Implicit velocity correction-based immersed boundary-lattice Boltzmann method and its applications, J. Comput. Phys. 228 (2009) 1963-1979.
  * [18] Y. Hu, H.Z. Yuan, S. Shu, X.D. Niu, M.J. Li, An improved momentum exchanged-based immersed boundary-lattice Boltzmann method by using an iterative technique, Comput. Math. Appl. 68 (3) (2014) 140-155.
  * [19] Y. Cai, S. Li, J. Lu, An improved immersed boundary-lattice Boltzmann method based on force correction technique, Int. J. Numer. Methods Fluids 87 (3) (2018) 109-133.
  * [20] D.R. Noble, J.R. Torczynski, A lattice-Boltzmann method for partially saturated computational cells, Int. J. Modern Phys. C. 9 (8) (1998) 1189-1201.
  * [21] B.K. Cook, D.R. Noble, J.R. Williams, A direct simulation method for particle-fluid systems, Eng. Comput. 21 (2-4) (2004) 151-168.
  * [22] T. Najuch, J. Sun, Analysis of two partially-saturated-cell methods for lattice Boltzmann simulation of granular suspension rheology, Comput. Fluids 189 (2019) 1-12.
  * [23] C. Tsigginos, J. Meng, X.J. Gu, D.R. Emerson, Lattice Boltzmann modeling of fluid-particle interaction based on a two-phase mixture representation, Phys. Rev. E 100 (2019) 063311.
  * [24] Y.H. Qian, D. d’Humières, P. Lallemand, Lattice BGK models for Navier-Stokes equation, Europhys. Lett. 17 (1992) 479-484.
  * [25] D.R.J. Owen, C.R. Leonardi, Y.T. Feng, An efficient framework for fluid-structure interaction using the lattice Boltzmann method and immersed moving boundaries, Int. J. Numer. Methods Eng. 87 (1-5) (2011) 66-95.
  * [26] J. Shen, Modeling and numerical approximation of two-phase incompressible flows by a phase-field approach, in: W. Bao, Q. Du (Eds.), Multiscale Modeling and Analysis for Materials Simulation, in: Lect. Note Ser., vol. 9, IMS, National University of Singapore, 2011, pp. 147-196.
  * [27] Z.L. Guo, C.G. Zheng, B.C. Shi, Non-equilibrium extrapolation method for velocity and pressure boundary conditions in the lattice Boltzmann method, Chin. Phys. 11 (4) (2002) 366.
  * [28] X.D. Niu, C. Shu, Y.T. Chew, Y. Peng, A momentum exchange-based immersed boundary-lattice Boltzmann method for simulating incompressible viscous flows, Phys. Lett. A 354 (3) (2006) 173-182.
  * [29] X. He, G. Doolen, Lattice Boltzmann method on curvilinear coordinates system: flow around a circular cylinder, J. Comput. Phys. 134 (2) (1997) 306-315.
  * [30] J. Happel, H. Brenner, _Low reynolds number hydrodynamics with special application to particulate media_ (Prentice-Hall, New York, 1965).
  * [31] L. Wang, Z. Guo, B. Shi, C. Zheng, Evaluation of three lattice Boltzmann models for particulate flows, Commun. Comput. Phys. 13 (4) (2013) 1151-1172.
  * [32] H. Li, X. Lu, H. Fang, Y. Qian, Force evaluations in lattice Boltzmann simulations with moving boundaries in two dimensions, Phys. Rev. E 70 (2) (2004) 026701.
  * [33] D. Wan, S. Turek, Direct numerical simulation of particulate flow via multigrid FEM techniques and the fictitious boundary method, Int. J. Numer. Methods Fluids 51 (2006) 531-66.
  * [34] S. Jafari, R. Yamamoto, M. Rahnama, Lattice-Boltzmann method combined with smoothed-profile method for particulate suspensions, Phys. Rev. E 83 (2011) 026702.
  * [35] J. Su, G. Chai, L. Wang, W. Cao, Z. Gu, C. Chen, X.Y. Xu, Pore-scale direct numerical simulation of particle transport in porous media, Chem. Eng. Sci. 199 (2019) 613-627.