# Robust Federated Learning via Over-The-Air Computation

Houssem Sifaou and Geoffrey Ye Li  Department of Electrical and Electronic
Engineering, Imperial College London, UK
emails<EMAIL_ADDRESS>
###### Abstract

This paper investigates the robustness of over-the-air federated learning to
Byzantine attacks. The simple averaging of the model updates via over-the-air
computation makes the learning task vulnerable to random or intended
modifications of the local model updates of some malicious clients. We propose
a robust transmission and aggregation framework to such attacks while
preserving the benefits of over-the-air computation for federated learning.
For the proposed robust federated learning, the participating clients are
randomly divided into groups and a transmission time slot is allocated to each
group. The parameter server aggregates the results of the different groups
using a robust aggregation technique and conveys the result to the clients for
another training round. We also analyze the convergence of the proposed
algorithm. Numerical simulations confirm the robustness of the proposed
approach to Byzantine attacks.

###### Index Terms:

Federated learning, Over-the-air computation, Byzantine attacks

##  I Introduction

With the rapid growth and unprecedented success of artificial intelligence
(AI) applications, a huge amount of data is transferred every day from
distributed clients to data centers for further processing. However, in many
applications, the data is sensitive, and sending it to data centers
constitutes a major privacy concern [1, 2, 3]. One promising solution is
federated learning (FL)[4, 5, 6]. With FL, multiple distributed clients train
a global machine learning model without sharing their local data. Local
computations are carried out at the different clients, and only the model
parameters updates are sent to the central server (CS). The latter aggregates
the local updates and forwards the result to the distributed clients for
another training round.

Although data privacy is preserved using FL, allowing the clients to perform
model parameter update opens the way for possible adversarial attacks. Some of
the distributed clients may send modified parameter updates with the intention
of misleading the learning process [7, 8]. In this context, Byzantine attacks
are a popular class, where certain clients aim to prevent the model from
converging or causing convergence to a falsified model. The Byzantine clients
may act independently or collectively. Unfortunately, even a single malicious
client in a distributed setup as in FL can seriously affect the end model [9].
Developing countermeasures to these attacks has gained an increasing interest
in recent years. Considering the FL setup without communication constraint,
several aggregation techniques have been proposed to robustify the stochastic
gradient descent (SGD) in the distributed setup. Stochastic algorithms
tolerating few Byzantine attacks has been developed by aggregating the
stochastic gradients updates using the median [10], geometric median[11, 12],
trimmed mean [13], and iterative filtering [14]. The Krum aggregation in [9]
selects the stochastic gradient having a minimum distance from a certain
number of nearest stochastic gradients. A robust stochastic aggregation (RSA)
has been developed in [15], which tolerates heterogeneous datasets and
Byzantine attacks. Other related works include leveraging redundancy of
gradients to enhance robustness [16, 17] and avoiding saddle points in non-
convex distributed learning in presence of Byzantine attacks [18]. Moreover,
the advantages of reducing the variance of the stochastic gradients to defend
against Byzantine attackers have been investigated in [19]. All these works
are applicable when the individual local updates are sent separately to the
parameter server. However, in the case of over-the-air FL (OTA-FL) [20, 21],
the local model updates are sent simultaneously over the analog wireless
channel. This makes the aforementioned robust aggregation techniques not
directly applicable in the case of OTA-FL.

In this work, we investigate the problem of Byzantine attacks in the case of
OTA-FL. Particularly, we propose a transmission and aggregation approach that
simultaneously exploits the benefits of over-the-air computation while being
robust to Byzantine attacks. By dividing the participating clients into
several groups at each global training round, assigning a transmission time
slot for each group, and aggregating the model updates of the different groups
using geometric median aggregation. Theoretical convergence analysis of the
proposed approach is conducted under some assumptions on the loss functions.
The analysis reveals that when the number of attacks is less than half of the
number of groups, the proposed algorithm converges at a linear rate to a
neighborhood of the optimal solution with an asymptotic learning error that
depends on the noise variance and the number of Byzantine attackers. Moreover,
as evidenced by numerical results, the proposed approach is robust to
Byzantine attacks compared with simple averaging OTA-FL.

The remainder of the paper is organized as follows. In the next section, we
introduce OTA-FL and the transmission model. In Section III, our proposed
robust approach is presented, while in Section IV the convergence of the
proposed algorithm is studied. Numerical results are provided in Section V and
concluding remarks are drawn in Section VI.

##  II System Model

We consider a FL system composed of a CS and $N$ clients. Client $n$ has its
local dataset $D_{n}=\\{({\bf
x}_{i}\in\mathbb{R}^{d},y_{i}\in\mathbb{R})\\}_{i=1}^{m_{n}}$ composed of
$m_{n}$ samples. All the clients collaboratively train a global model by
communicating with the CS. More precisely, they seek for the optimal parameter
vector ${\bf w}^{\star}\in\mathbb{R}^{p}$ that minimizes a global loss
function $f({\bf w})$ given by

$f({\bf w})=\frac{1}{N}\sum_{n=1}^{N}f_{n}({\bf w}),$ (1)

where $f_{n}({\bf w})$ is the local loss function at client $n$ defined as

$f_{n}({\bf w})=\frac{1}{m_{n}}\sum_{j=1}^{m_{n}}\ell({\bf w};{\bf
x}_{j},y_{j}).$ (2)

Usually, the CS communicates with the clients using wireless channels. We
consider in this work the case of analog OTA-FL [21], which will be introduced
hereafter.

###  II-A Analog over-the-air FL

At each global training round $t$, the CS sends the model parameter vector,
${\bf w}_{t}$, to the clients. It is usually assumed that the downlink
communication is perfect due to the high power available at the CS. Thus, each
client receives the global model without distortions. Then, client $n$ sets
its local model as ${\bf w}_{t,0}^{n}={\bf w}_{t}$ and runs its local SGD for
$H$ iterations based on its local dataset

${\bf w}_{t,i+1}^{n}={\bf w}_{t,i}^{n}-\eta_{t}f_{n,i_{n}^{t}}^{\prime}({\bf
w}_{t,i}^{n}),\ \ {\rm for}\ \ i=0,1,\cdots,H-1,$ (3)

where $\eta_{t}$ is the SGD step size at round $t$ and
$f_{n,i_{n}^{t}}^{\prime}({\bf w}_{t,i}^{n})$ denotes the stochastic gradient
computed using a sample with index $i_{n}^{t}$ chosen uniformly at random from
the local dataset of client $n$. In practice, After $H$ iterations, the
clients convey their model updates given by

${\bf m}_{t}^{n}={\bf w}_{t,H}^{n}-{\bf w}_{t},$ (4)

simultaneously to the CS via analog OTA-FL. The model updates should be
precoded in order to mitigate the effect of channel fading. Let
$\tilde{h}_{n,t}=h_{t,n}e^{j\phi_{t}^{n}}$ be the a block fading channel
corresponding to user $n$ at the transmission time of global round $t$, where
$h_{t,n}>0$ and $\phi_{t}^{n}=[-\pi,\pi]$ are its module and phase
respectively. As in [20, 22, 21, 23, 24], we assume that perfect channel state
information (CSI) is available at the clients and the CS. The imperfect CSI
case is left for future investigation. Moreover, since the power budget at the
clients is limited, the transmitted signal should satisfy the following
average power constraint

$\mathbb{E}\left[\|{\bf x}_{n}\|^{2}\right]\leq P.$ (5)

In practice, weak channels might cause a high amplification of transmit power,
possibly violating the transmission power constraint (5). To overcome this
issue, a threshold $h_{\min}$ can be set, and clients with channel fading
coefficients less than $h_{\min}$ in magnitude will not transmit in that
training round.

We adopt in this paper the precoding scheme proposed in [23], where every
client $n$ precodes its model update ${\bf m}_{t}^{n}$ as

${\bf x}_{t,n}=\begin{cases}\rho_{t}\frac{h_{\rm
min}}{h_{t,n}}e^{-j\phi_{t}^{n}}{\bf m}_{t}^{n},\ \ \ {\rm if}\ \
h_{t,n}>h_{\min}\\\ 0,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {\rm if}\ \
h_{t,n}\leq h_{\min}\end{cases}$ (6)

where $\rho_{t}$ is a factor to satisfy the average power constraint at client
$n$. As proposed in [23], $\rho_{t}$ can be set as follows

$\rho_{t}=\sqrt{\frac{P}{\max_{n}\mathbb{E}\|{\bf m}_{t}^{n}\|^{2}}}.$ (7)

Note that in practice, the clients do not have access the updates of each
other in order to compute $\rho_{t}$ at each training round. A possible way to
deal with this issue is that the CS can estimate this parameter offline using
a small dataset and then forward it to the clients so they can use it at every
training round. Another solution is to obtain an upper bound for $\rho_{t}$
and use it at every iteration [23]. The received signal at the CS is

${\bf y}_{t}=\sum_{n\in\mathcal{K}_{t}}\rho_{t}h_{\min}{\bf
m}^{n}_{t}+\tilde{\bf z}_{t},$ (8)

where $\mathcal{K}_{t}$ is the set of clients indices with channel fading
verifying $h_{n,t}>h_{\min}$ and $\tilde{\bf
z}_{t}\sim\mathcal{N}(\boldsymbol{0},\sigma^{2}{\bf I}_{p})$ stands for
additive noise. In order to update the global model, the CS sets

${\bf w}_{t+1}=\frac{{\bf y}_{t}}{|\mathcal{K}_{t}|\rho_{t}h_{\min}}+{\bf
w}_{t},$ (9)

where $|\mathcal{K}_{t}|$ is the cardinality of the set $\mathcal{K}_{t}$. The
global model update in (9) can be also written as

${\bf w}_{t+1}=\frac{1}{|\mathcal{K}_{t}|}\sum_{n=1}^{|\mathcal{K}_{t}|}{\bf
w}^{n}_{t,H}+{\bf z}_{t},$ (10)

where ${\bf z}_{t}\triangleq\frac{\tilde{\bf
z}_{t}}{|\mathcal{K}_{t}|\rho_{t}h_{\min}}\sim\mathcal{N}(\boldsymbol{0},\frac{\sigma^{2}}{|\mathcal{K}_{t}|^{2}\rho_{t}^{2}h_{\min}^{2}}{\bf
I}_{p})$.

###  II-B Byzantine attacks

We assume that $B<N$ clients are malicious; sending arbitrary or modified
parameter vector updates aiming at affecting the convergence of the global
model or forcing it to converge to some particular solution. This type of
attacks is known as Byzantine attacks. There are many works that proposed
solutions to deal with this type of attacks in federated learning [9, 10, 11,
12, 13, 14, 15, 16, 17, 18, 19]. However, these works considered the case of
wired FL or the case that each individual update of the client is sent
separately to the CS. In this work, we consider the effect of Byzantine
attacks in the context of analog OTA-FL and we propose a practical approach to
deal with such attacks.

##  III Proposed approach

In order to reduce the effect of the malicious model updates, we propose the
following approach. At each global training round $t$, the CS divides
uniformly at random the $N$ clients into $G=N/m$ groups where each group is
composed of $m$ clients. Each group will be allocated a time slot for
transmission of their model updates. Precisely, the clients of group $g$ will
transmit simultaneously their updates over-the-air. This allows the CS to
obtain $G$ model updates. Then, with a robust aggregation technique, the
different model updates of the groups will be aggregated to update the global
model. This approach will be robust to Byzantine attacks.

At global iteration $t$, the global model, ${\bf w}_{t}$, is forwarded to all
the clients. The clients in group $g$ perform $H_{g}$ steps of SGD using their
local datasets as in $\eqref{SGD}$, and compute their model updates

${\bf m}_{t}^{n}={\bf w}_{t,H_{g}}^{n}-{\bf w}_{t},\ \ {\rm for}\ \
n\in\mathcal{G}_{g,t},$ (11)

where $\mathcal{G}_{g,t}$ is the set of user indices belonging to group $g$ at
the global training round $t$. Note that since the clients in different groups
are not sending at the same time slot, we can let the number of the SGD steps
varying among groups. In other terms, the clients of a group with later
transmission time can perform more SGD steps than those in the current
transmitting group. However, for simplicity we assume in the sequel that all
clients perform the same number of SGD steps $H$ regardless of their
transmission time. The clients in group $g$ compute their precoded signal as

${\bf x}_{t,n}=\begin{cases}\rho_{t}\frac{h_{\rm
min}}{h_{t,n}}e^{-j\phi_{t}^{n}}{\bf m}_{t}^{n},\ \ \ {\rm if}\ \
h_{t,n}>h_{\min}\\\ 0,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {\rm if}\ \
h_{t,n}\leq h_{\min}\end{cases}$ (12)

and send their updates simultaneously during their allocated transmission time
slot. At the CS, the received signal vector corresponding to group $g$ can be
expressed as

${\bf y}_{t,g}=\sum_{n\in\mathcal{K}_{t,g}}\rho_{t}h_{\min}{\bf
m}^{n}_{t}+\tilde{\bf z}_{t,g},$ (13)

where $\mathcal{K}_{t,g}$ is the set of clients indices of group $g$ with
channels such that $h_{n,t}>h_{\min}$ and $\tilde{\bf
z}_{t,g}\sim\mathcal{N}(\boldsymbol{0},\sigma^{2}{\bf I}_{p})$ is the additive
noise. The CS estimates the model update of the group $g$ as

${\bf u}_{g}^{t}=\frac{{\bf y}_{t,g}}{\rho_{t}h_{\min}|\mathcal{K}_{t,g}|}.$
(14)

After all the group updates are collected, the CS disposes of $G$ vector
updates ${\bf u}_{1}^{t},\cdots,{\bf u}_{G}^{t}$ and can aggregate these
updates to obtain the global model. For instance, one of the most efficient
aggregation techniques that can be used is the geometric median [25]. Other
aggregation techniques can be used such as the Krum aggregation rule proposed
in [9]. In this work, we focus on geometric median aggregation. The global
model is updated as

${\bf w}_{t+1}={\rm geomed}({\bf u}_{1}^{t},\cdots,{\bf u}_{G}^{t})+{\bf
w}_{t},$ (15)

where ${\rm geomed}(.)$ stands for the geometric median aggregation defined as

${\rm geomed}(\\{{\bf u}_{i}\\}_{i\in\mathcal{K}})={\rm arg}\min_{\bf
z}\sum_{i\in\mathcal{K}}\|{\bf z}-{\bf u}_{i}\|.$

The geometric median aggregation has been proposed as an efficient solution to
Byzantine attacks when the individual updates of the clients are sent
separately to the CS. In fact, it approximates well the mean of the honest
clients weight updates when $B<N/2$ [19]. In our case, we will use it to
aggregate the group updates. The number of Byzantine workers should satisfy
$B<G/2$ in order for the geometric median to well approximate the mean of the
group updates composed by only honest clients.

To compute the geometric median, the Weiszfeld algorithm can be used [26]. To
avoid numerical instabilities, a smoothed version of the Weiszfeld algorithm
can be implemented in practice [27], which computes a smoothed geometric
median defined as ${\rm geomed}_{\epsilon}(\\{{\bf
u}_{i}\\}_{i\in\mathcal{K}})={\rm arg}\min_{\bf z}\sum_{i\in\mathcal{K}}\|{\bf
z}-{\bf u}_{i}\|_{\epsilon},$ where

$\|{\bf x}\|_{\epsilon}=\begin{cases}\frac{1}{2\epsilon}\|{\bf
x}\|^{2}+\frac{\epsilon}{2}\ \ \ {\rm if}\ \ \|{\bf x}\|\leq\epsilon\\\ \|{\bf
x}\|\ \ \ \ \ \ \ \ \ \ \ \ \ {\rm if}\ \ \|{\bf x}\|>\epsilon,\end{cases}$

where $\epsilon>0$ is a smoothing parameter. The steps of the proposed
approach are summarised in Algorithm 1, where $LocalComp({\bf
w}_{t},H,b,\eta,D_{n})$ consists of $H$ steps of batch-SGD using local dataset
$D_{n}$ with learning rate $\eta$ and initial parameter vector ${\bf w}_{t}$
as described in (3).

Input: Initial model ${\bf w}_{0}$

for  _$t=0,1,2,\cdots$_ do

The CS forwards ${\bf w}_{t}$ to the clients;

for  _each client $n$_ do

${\bf m}_{t}^{n}\leftarrow LocalComp({\bf w}_{t},H,b,\eta,D_{n})$

end for

for  _$g=1,2,\cdots,G$_ do

for  _each client $n$ in group $g$ $(n\in\mathcal{G}_{g,t})$_ do

client $n$ transmit its model update ${\bf x}_{t,n}$ precoded via (12) during
transmission time slot $T_{g}$

end for

The CS receives ${\bf y}_{g}$ of group $g$ and computes ${\bf u}_{t}^{g}$ as
in (14)

end for

The CS aggregates the group updates using (15) to obtain the new global
parameter vector ${\bf w}_{t+1}$

end for

Algorithm 1

##  IV Convergence Analysis

In this section, we analyze the convergence of the proposed approach. Several
assumptions on the loss functions are needed.

######  Assumption 1.

  * (i)

Strong convexity: The objective function $f$ is $\mu$-strongly convex, that
is, for all ${\bf x},{\bf y}\in\mathbb{R}^{p}$

$f({\bf x})\geq f({\bf y})+\langle f^{\prime}({\bf y}),{\bf x}-{\bf
y}\rangle+\frac{\mu}{2}\|{\bf x}-{\bf y}\|^{2}.$

for some $\mu>0$.

  * (ii)

Lipschitz continuity of gradients: The objective function $f$ has $L$-Lipchitz
continuous gradients, that is, for all ${\bf x},{\bf y}\in\mathbb{R}^{p}$

$\|f^{\prime}({\bf x})-f^{\prime}({\bf y})\|\leq L\|{\bf x}-{\bf y}\|.$

for some $L>0$.

  * (iii)

Bounded outer variation: For each honest client $n$, the variations of its
aggregated gradients with respect to the over-all gradient is bounded as

$\mathbb{E}\|f^{\prime}_{n}({\bf x})-f^{\prime}({\bf x})\|^{2}\leq\delta^{2},\
{\rm for\ all}\ {\bf x}\in\mathbb{R}^{p}.$

  * (iv)

Bounded inner variation: for each honest client $n$, the variation of its
stochastic gradients with respect to its aggregated gradients is bounded as

$\mathbb{E}\|f^{\prime}_{n,i_{n}^{t}}({\bf x})-f_{n}^{\prime}({\bf
x})\|^{2}\leq\kappa^{2},\ {\rm for\ all}\ {\bf x}\in\mathbb{R}^{p}.$

  * (v)

Bounded stochastic gradients: For each honest client $n$, stochastic gradient
$f^{\prime}_{n,i_{n}^{t}}({\bf x})$ satisfies

$\mathbb{E}\|f^{\prime}_{n,i_{n}^{t}}({\bf x})\|^{2}\leq K^{2},\ {\rm for\
all}\ {\bf x}\in\mathbb{R}^{p},$

for some fixed $K^{2}>0$.

Items (i) and (ii) in Assumption 1 are common in convex analysis. Items (iii)
and (iv) are needed to bound the inner and outer variations of the stochastic
gradients and the gradients of the honest clients, respectively. These
assumptions are adopted in most of the existing works considering distributed
SGD in presence of Byzantine attacks [19, 28]. The convergence of the proposed
approach is presented in the following theorem. For simplicity, we assume in
this section that the learning rate is constant and that the clients perform
one SGD step at each global iteration, that is, $H=1$ and $\eta_{t}=\eta$ for
all $t$.

######  Theorem 1.

Under Assumption 1, when the number of Byzantine attackers satisfies
$B<\frac{G}{2}$ and the step size $\eta$ verifies $\eta<\frac{\mu}{2L^{2}}$,
then

$\mathbb{E}\|{\bf w}_{t}-{\bf w}^{*}\|^{2}\leq(1-\eta\mu)^{t}B+A,$

where

$\displaystyle B$ $\displaystyle=\|{\bf w}_{0}-{\bf w}^{*}\|^{2}-A,$
$\displaystyle A$
$\displaystyle=\frac{2}{\mu^{2}}C_{\alpha}^{2}\left(\delta^{2}+\kappa^{2}+\frac{p\sigma^{2}}{Ph_{min}^{2}}K^{2}\right),$

with $C_{\alpha}=\frac{2-2\alpha}{1-2\alpha}$ and $\alpha=\frac{B}{G}$.

###### Proof:

See Appendix. ∎

Theorem 1 states that the proposed approach converges at a linear rate to a
neighborhood of ${\bf w}^{*}$. The asymptotic learning error, $A$, depends on
the number of Byzantine attackers through $C_{\alpha}$ and on the noise
variance. When the number of Byzantine attackers increases, $C_{\alpha}$
increases, which yields a higher asymptotic error.

##  V Numerical Results

In order to evaluate the performance of the proposed approach, we provide in
this section simulation on real dataset. We consider the MNIST dataset
composed of $28\times 28$ images of handwritten digits corresponding to 10
classes (digits 0 to 9). The dataset set contains $60,000$ training samples
and 10,000 testing samples. We divided the training samples equally at random
to $K=100$ clients, that is, each client has a local dataset composed of $600$
samples. We used the multi-class logistic regression model. At every global
training round, each client performs $H=12$ steps of local batch-SGD where at
each step a minibatch of size $b=50$ is used. The leaning rate $\eta_{t}$ is
fixed to $0.01$ in all local SGD steps.The number of groups is fixed to
$G=20$, while the noise variance and the power constraint are set to be
$\sigma^{2}=10^{-2}$ and $P=1$, respectively. The smoothing parameter of the
Weiszfeld algorithm used is $\epsilon=10^{-4}$. In all simulations, the
transmission threshold $h_{\min}$ and the scaling factor $\rho_{t}$ are fixed
to 0.1 and 8 respectively.

In the first experiment, we consider Gaussian attacks where each Byzantine
client sends a Gaussian vector of mean ${\bf w}_{t}$ and variance $30$ instead
of its actual model update. We compare our proposed approach and the simple
averaging COTAF [23] described in Section II. From Fig. 1, simple averaging is
not robust to Byzantine attacks. On the other hand, our proposed approach
guarantees fast convergence.

In the second experiment, we consider class flip attacks where the Byzantine
clients change the labels of their local datasets as $y=9-y$. Fig. 2
demonstrates the effect of increasing the number of Byzantine workers. From
this figure, the proposed approach is robust to Byzantine attacks even though
the test accuracy relatively decreases with the number of the Byzantine
clients. This is expected since the training, in presence of malicious
clients, is done over a smaller sample size as only the honest workers
contribute in the learning process. This was also predicted by the convergence
analysis conducted in the previous section, where it has been shown that the
asymptotic learning error increases with the number of Byzantine attacks.

$0$$100$$200$$300$$400$$500$$10^{-1}$$10^{0}$$10^{1}$$10^{2}$IterationTest
LossProposed ApproachSimple averaging (a) Test Loss

$0$$100$$200$$300$$400$$500$$0$$0.2$$0.4$$0.6$$0.8$$1$IterationTest
AccuracyProposed ApproachSimple averaging (b) Test Accuracy

Figure 1: Test loss and test accuracy vs. the number of iterations when $B=9$
Byzantine clients applying Gaussian attacks.

$0$$100$$200$$300$$400$$500$$10^{-1}$$10^{0}$$10^{1}$IterationTest Loss$B=10$
$B=5$ $B=0$ (a) Test Loss

$0$$100$$200$$300$$400$$500$$0$$0.2$$0.4$$0.6$$0.8$$1$IterationTest
Accuracy$B=10$ $B=5$ $B=0$ (b) Test Accuracy

Figure 2: Test loss and test accuracy vs. the number of iterations for
different numbers of Byzantine clients applying class-flip attacks.

##  VI Conclusion

In this paper, we have proposed a novel approach to account for Byzantine
attacks in over-the-air FL. By dividing the distributed clients into groups,
the parameter server is able to reduce the effect of Byzantine attacks by
aggregating the group parameter updates. The convergence of the proposed
algorithm has been studied analytically under assumption of convex loss
functions. The simulation results show the robustness of the approach to
different Byzantine attacks. This work can be extended by studying other
aggregation techniques. It is also important to study the effect of the
imperfect channel knowledge on the proposed algorithm.

## Appendix

We first sate the following lemma which will be used later in the proof.

######  Lemma 2.

[19, Lemma 2] Let $\mathcal{V}$ be a subset of random vectors distributed in a
normed vector space. If $\mathcal{V}^{\prime}\subset\mathcal{V}$ such that
$|\mathcal{V}^{\prime}|<\frac{|\mathcal{V}|}{2}$, then it holds that

$\mathbb{E}\|\underset{{\bf v}\in\mathcal{V}}{{\rm geomed}}({\bf v})\|^{2}\leq
C_{\alpha}^{2}\frac{\sum_{{\bf v}\notin\mathcal{V}^{\prime}}\|{\bf
v}\|^{2}}{|\mathcal{V}|-|\mathcal{V}^{\prime}|},$

where $C_{\alpha}=\frac{2-2\alpha}{1-2\alpha}$ and
$\alpha=\frac{|\mathcal{V}^{\prime}|}{|\mathcal{V}|}$.

Define $\delta_{t}\triangleq\mathbb{E}\|{\bf w}_{t}-{\bf w}^{*}\|^{2}$. To
prove Theorem 1, we start by finding an upper bound for $\delta_{t+1}$,

$\displaystyle\delta_{t+1}$ $\displaystyle=\mathbb{E}\|{\bf w}_{t+1}-{\bf
w}^{*}\|^{2}$ $\displaystyle=\mathbb{E}\|{\bf w}_{t}-\eta f^{\prime}({\bf
w}_{t})-{\bf w}^{*}+{\bf w}_{t+1}-{\bf w}_{t}+\eta f^{\prime}({\bf
w}_{t})\|^{2}$ $\displaystyle\leq\frac{1}{1-\gamma}\mathbb{E}\|{\bf
w}_{t}-\eta f^{\prime}({\bf w}_{t})-{\bf w}^{*}\|^{2}$ $\displaystyle\ \ \ \
+\frac{1}{\gamma}\mathbb{E}\|{\bf w}_{t+1}-{\bf w}_{t}+\eta f^{\prime}({\bf
w}_{t})\|^{2},$

for any $0<\gamma<1$, where we have used the fact that $\|{\bf x}+{\bf
y}\|^{2}\leq\frac{1}{1-\gamma}\|{\bf x}\|^{2}+\frac{1}{\gamma}\|{\bf y}\|^{2}$
for any $0<\gamma<1$. Since $f^{\prime}({\bf w}^{*})=0$, we can write

$\displaystyle\|{\bf w}_{t}-\eta f^{\prime}({\bf w}_{t})-{\bf w}^{*}\|^{2}$
$\displaystyle=\|{\bf w}_{t}-\eta(f^{\prime}({\bf w}_{t})-f^{\prime}({\bf
w}^{*}))-{\bf w}^{*}\|^{2}$ $\displaystyle=\|{\bf w}_{t}-{\bf
w}^{*}\|^{2}-2\eta\langle f^{\prime}({\bf w}_{t})-f^{\prime}({\bf w}^{*}),{\bf
w}_{t}-{\bf w}^{*}\rangle$ $\displaystyle\ \ \ \ \ +\eta^{2}\|f^{\prime}({\bf
w}_{t})-f^{\prime}({\bf w}^{*})\|^{2}$ $\displaystyle\overset{(a)}{\leq}\|{\bf
w}_{t}-{\bf w}^{*}\|^{2}-2\eta\mu\|{\bf w}_{t}-{\bf
w}^{*}\|^{2}+\eta^{2}L^{2}\|{\bf w}_{t}-{\bf w}^{*}\|^{2}$

where $(a)$ follows from items $(i)$ and $(ii)$ of Assumption 1. Thus,

$\displaystyle\delta_{t+1}\leq\frac{1-2\eta\mu+\eta^{2}L^{2}}{1-\gamma}\delta_{t}+\frac{1}{\gamma}\mathbb{E}\|{\bf
w}_{t+1}-{\bf w}_{t}+\eta f^{\prime}({\bf w}_{t})\|^{2}.$

For $\eta<\frac{2}{\mu}$, we can take $\gamma=\frac{\eta\mu}{2}$. Assuming
further that $\eta\leq\frac{\mu}{2L^{2}}$, it holds that
$\frac{1-2\eta\mu+\eta^{2}L^{2}}{1-\gamma}\leq 1-\eta\mu$. Hence, for
$\eta<\min(\frac{2}{\mu},\frac{\mu}{2L^{2}})=\frac{\mu}{2L^{2}}$,

$\displaystyle\delta_{t+1}$
$\displaystyle\leq(1-\eta\mu)\delta_{t}+\frac{2}{\eta\mu}\mathbb{E}\|{\bf
w}_{t+1}-{\bf w}_{t}+\eta f^{\prime}({\bf w}_{t})\|^{2}.$ (16)

We treat now the second term of the right-hand side of (16). From the update
rule (15), it follows that ${\bf w}_{t+1}-{\bf w}_{t}={\rm geomed}(\\{{\bf
u}_{g}^{t}\\}_{g=1}^{G})$ where ${\bf u}_{g}^{t}$, for the case $H=1$, is
given by

$\displaystyle{\bf u}_{g}^{t}$
$\displaystyle=-\frac{1}{|\mathcal{K}_{t,g}|}\sum_{n\in\mathcal{K}_{t,g}}\eta
f_{n,i_{n}^{t}}^{\prime}({\bf w}_{t})+{\bf z}_{t,g},$

Define $\mathcal{B}_{t}$ as the set of groups containing at least one
Byzantine attacker at global iteration $t$. Applying Lemma 2 yields

$\displaystyle\mathbb{E}\|{\bf w}_{t+1}-{\bf w}_{t}+\eta f^{\prime}({\bf
w}_{t})\|^{2}=\mathbb{E}\|{\rm geomed}(\\{{\bf u}_{g}^{t}\\}_{g=1}^{G})+\eta
f^{\prime}({\bf w}_{t})\|^{2}$ $\displaystyle=\mathbb{E}\|{\rm geomed}(\\{{\bf
u}_{g}^{t}+\eta f^{\prime}({\bf w}_{t})\\}_{g=1}^{G})\|^{2}$
$\displaystyle\leq
C_{\alpha}^{2}\frac{\sum_{g\notin\mathcal{B}_{t}}\mathbb{E}\left\|{\bf
u}_{g}^{t}+\eta f^{\prime}({\bf w}_{t})\right\|^{2}}{G-|\mathcal{B}_{t}|},$

Replacing ${\bf u}_{g}^{t}$ by its expression, it follows that

$\displaystyle\mathbb{E}\|{\bf w}_{t+1}-{\bf w}_{t}+\eta f^{\prime}({\bf
w}_{t})\|^{2}$ $\displaystyle\leq
C_{\alpha}^{2}\frac{\sum_{g\notin\mathcal{B}_{t}}\mathbb{E}\left\|-\frac{\eta}{|\mathcal{K}_{t,g}|}\sum_{n\in\mathcal{K}_{t,g}}f^{\prime}_{i^{n}}({\bf
w}_{t})+{\bf z}_{t,g}+\eta f^{\prime}({\bf
w}_{t})\right\|^{2}}{G-|\mathcal{B}_{t}|}$ $\displaystyle\leq
C_{\alpha}^{2}\eta^{2}\frac{\sum_{g\notin\mathcal{B}_{t}}\frac{1}{|\mathcal{K}_{t,g}|}\sum_{n\in\mathcal{K}_{t,g}}\mathbb{E}\left\|f^{\prime}_{n,i_{n}^{t}}({\bf
w}_{t})-f^{\prime}_{n}({\bf w}_{t})\right\|^{2}}{G-|\mathcal{B}_{t}|}$
$\displaystyle+C_{\alpha}^{2}\eta^{2}\frac{\sum_{g\notin\mathcal{B}_{t}}\frac{1}{|\mathcal{K}_{t,g}|}\sum_{n\in\mathcal{K}_{t,g}}\mathbb{E}\left\|f^{\prime}_{n}({\bf
w}_{t})-f^{\prime}({\bf w}_{t})\right\|^{2}}{G-|\mathcal{B}_{t}|}$
$\displaystyle+C_{\alpha}^{2}\frac{\sum_{g\notin\mathcal{B}_{t}}\mathbb{E}\|{\bf
z}_{t,g}\|^{2}}{G-|\mathcal{B}_{t}|},$

where the last result follows from the fact that $\left\|\sum_{i=1}^{k}{\bf
v}_{i}\right\|^{2}\leq k\sum_{i=1}^{k}\|{\bf v}_{i}\|^{2}$ for any sequence of
vectors $\\{{\bf v}_{i}\\}_{i=1}^{k}$ and the independence between the noise
vectors and the stochastic gradients. $\mathbb{E}\|{\bf z}_{t,g}\|^{2}$ can be
bounded as follows

$\displaystyle\mathbb{E}\|{\bf
z}_{t,g}\|^{2}=\frac{p\sigma^{2}}{h_{min}^{2}|\mathcal{K}_{t,g}|^{2}\rho_{t}^{2}}=\frac{p\sigma^{2}}{h_{min}^{2}|\mathcal{K}_{t,g}|^{2}}\frac{\max_{n}\mathbb{E}\|{\bf
m}_{n}^{t}\|^{2}}{P}$
$\displaystyle\leq\frac{p\sigma^{2}}{Ph_{min}^{2}|\mathcal{K}_{t,g}|^{2}}\eta^{2}K^{2}\leq\frac{p\sigma^{2}}{Ph_{min}^{2}}\eta^{2}K^{2}.$
(17)

Combining (17) and items (iii) and (iv) in Assumption 1, it holds that

$\displaystyle\mathbb{E}\|{\bf w}_{t+1}-{\bf w}_{t}+\eta f^{\prime}({\bf
w}_{t})\|^{2}\\!\leq
C_{\alpha}^{2}\eta^{2}\\!\left(\delta^{2}+\kappa^{2}+\frac{p\sigma^{2}}{Ph_{min}^{2}}K^{2}\right)\\!.$
(18)

Combining (16) and (18) yields for $\eta<\frac{\mu}{2L^{2}}$

$\delta_{t+1}\leq(1-\eta\mu)\delta_{t}+\eta\mu A,$

where
$A\triangleq\frac{2}{\mu^{2}}C_{\alpha}^{2}\left(\delta^{2}+\kappa^{2}+\frac{p\sigma^{2}}{Ph_{min}^{2}}K^{2}\right)$.
Thus,

$\displaystyle\delta_{t+1}$
$\displaystyle\leq(1-\eta\mu)^{t+1}\delta_{0}+\eta\mu
A\sum_{i=0}^{t}(1-\eta\mu)^{i}$ (19)
$\displaystyle=(1-\eta\mu)^{t+1}(\delta_{0}-A)+A,$ (20)

which completes the proof.

## References

  * [1] R. Agrawal and R. Srikant, “Privacy-preserving data mining,” _SIGMOD Rec._ , vol. 29, no. 2, pp. 439–450, 2000.
  * [2] J. C. Duchi, M. I. Jordan, and M. J. Wainwright, “Local privacy and minimax bounds: sharp rates for probability estimation,” in _Proc. NeurIPS_ , 2013, pp. 1529–1537.
  * [3] L. Zhou, K.-H. Yeh, G. Hancke, Z. Liu, and C. Su, “Security and privacy for the industrial internet of things: An overview of approaches to safeguarding endpoints,” _IEEE Signal Processing Mag._ , vol. 35, no. 5, pp. 76–87, 2018\.
  * [4] B. McMahan, E. Moore, D. Ramage, S. Hampson, and B. A. y. Arcas, “Communication-efficient learning of deep networks from decentralized data,” in _Proc. AISTATS_ , vol. 54, 2017, pp. 1273–1282.
  * [5] P. Kairouz, H. B. McMahan, B. Avent, A. Bellet, M. Bennis, and A. N. B. et al., “Advances and open problems in federated learning,” _CoRR_ , vol. abs/1912.04977, 2019. [Online]. Available: http://arxiv.org/abs/1912.04977
  * [6] Z. Qin, G. Ye Li, and H. Ye, “Federated learning and wireless communications,” _IEEE Wireless Commun._ , pp. 1–7, 2021.
  * [7] A. Vempaty, L. Tong, and P. K. Varshney, “Distributed inference with byzantine data: State-of-the-art review on data falsification attacks,” _IEEE Signal Processing Mag._ , vol. 30, no. 5, pp. 65–75, 2013.
  * [8] Y. Chen, S. Kar, and J. M. Moura, “The internet of things: Secure distributed inference,” _IEEE Signal Processing Mag._ , vol. 35, no. 5, pp. 64–75, 2018\.
  * [9] P. Blanchard, E. M. El Mhamdi, R. Guerraoui, and J. Stainer, “Machine learning with adversaries: Byzantine tolerant gradient descent,” in _Proc. NeurIPS_ , 2017, pp. 118–128.
  * [10] C. Xie, O. Koyejo, and I. Gupta, “Generalized byzantine-tolerant sgd,” _ArXiv_ , vol. abs/1802.10116, 2018.
  * [11] S. Minsker, “Geometric median and robust estimation in banach spaces,” _Bernoulli_ , vol. 21, no. 4, pp. 2308–2335, 2015.
  * [12] Y. Chen, L. Su, and J. Xu, “Distributed statistical machine learning in adversarial settings: Byzantine gradient descent,” _Proc. ACM Meas. Anal. Comput. Syst._ , vol. 1, no. 2, 2017.
  * [13] D. Yin, Y. Chen, R. Kannan, and P. Bartlett, “Byzantine-robust distributed learning: Towards optimal statistical rates,” in _Proc. ICML_ , 2018, pp. 5650–5659.
  * [14] L. Su and J. Xu, “Securing distributed gradient descent in high dimensional statistical learning,” _Proc. ACM Meas. Anal. Comput. Syst._ , vol. 3, no. 1, 2019.
  * [15] L. Li, W. Xu, T. Chen, G. B. Giannakis, and Q. Ling, “Rsa: Byzantine-robust stochastic aggregation methods for distributed learning from heterogeneous datasets,” in _Proc. AAAI_ , vol. 33, no. 01, 2019, pp. 1544–1551.
  * [16] L. Chen, H. Wang, Z. Charles, and D. Papailiopoulos, “Draco: Byzantine-resilient distributed training via redundant gradients,” in _Proc. ICML_ , 2018, pp. 903–912.
  * [17] S. Rajput, H. Wang, Z. B. Charles, and D. Papailiopoulos, “Detox: A redundancy-based framework for faster and more robust gradient aggregation,” in _Proc. NeurIPS_ , 2019.
  * [18] D. Yin, Y. Chen, R. Kannan, and P. Bartlett, “Defending against saddle point attack in byzantine-robust distributed learning,” in _Proc. ICML_ , 2019, pp. 7074–7084.
  * [19] Z. Wu, Q. Ling, T. Chen, and G. B. Giannakis, “Federated variance-reduced stochastic gradient descent with robustness to byzantine attacks,” _IEEE Trans. Signal Process._ , vol. 68, pp. 4583–4596, 2020.
  * [20] K. Yang, T. Jiang, Y. Shi, and Z. Ding, “Federated learning via over-the-air computation,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 3, pp. 2022–2035, 2020.
  * [21] M. Mohammadi Amiri and D. Gündüz, “Machine learning at the wireless edge: Distributed stochastic gradient descent over-the-air,” _IEEE Trans. Signal Process._ , vol. 68, pp. 2155–2169, 2020.
  * [22] T. Sery and K. Cohen, “On analog gradient descent learning over multiple access fading channels,” _IEEE Transactions on Signal Processing_ , vol. 68, pp. 2897–2911, 2020.
  * [23] T. Sery, N. Shlezinger, K. Cohen, and Y. C. Eldar, “Over-the-air federated learning from heterogeneous data,” _CoRR_ , vol. abs/2009.12787, 2020. [Online]. Available: https://arxiv.org/abs/2009.12787
  * [24] D. Liu and O. Simeone, “Privacy for free: Wireless federated learning via uncoded transmission with adaptive power control,” _IEEE J. Sel. Areas Commun_ , vol. 39, no. 1, pp. 170–185, 2021.
  * [25] A. N. Bhagoji, S. Chakraborty, P. Mittal, and S. Calo, “Analyzing federated learning through an adversarial lens,” in _Proc. ICML_ , 2019, pp. 634–643.
  * [26] E. Weiszfeld, “Sur le point pour lequel la somme des distances de n points donnés est minimum,” _Tohoku Mathematical Journal, First Series_ , vol. 43, pp. 355–386, 1937.
  * [27] K. Pillutla, S. M. Kakade, and Z. Harchaoui, “Robust aggregation for federated learning,” 2019.
  * [28] H. Tang, X. Lian, M. Yan, C. Zhang, and J. Liu, “${D}^{2}$: Decentralized training over decentralized data,” in _Proc. ICML_ , J. Dy and A. Krause, Eds., vol. 80, 2018, pp. 4848–4856.


# Volume filtered FEM-DEM framework for simulating particle-laden flows in
complex geometries

Abhilash Reddy Malipeddi Corresponding author. e-mail<EMAIL_ADDRESS>Department of Mechanical Engineering, University of Michigan, Ann Arbor, MI
48109, USA. C. Alberto Figueroa  Department of Biomedical Engineering,
University of Michigan, Ann Arbor, MI 48109, USA.  Department of Surgery,
University of Michigan, Ann Arbor, MI 48109, USA.  Jesse Capecelatro
Department of Mechanical Engineering, University of Michigan, Ann Arbor, MI
48109, USA. Department of Aerospace Engineering, University of Michigan, Ann
Arbor, MI 48109, USA.

###### Abstract

We present a computational framework for modeling large-scale particle-laden
flows in complex domains with the goal of enabling simulations in medical-
image derived patient specific geometries. The framework is based on a volume-
filtered Eulerian-Lagrangian strategy that uses a finite element method (FEM)
to solve for the fluid phase coupled with a discrete element method (DEM) for
the particle phase, with varying levels of coupling between the phases. The
fluid phase is solved on a three-dimensional unstructured grid using a
stabilized FEM. The particle phase is modeled as rigid spheres and their
motion is calculated according to Newton’s second law for translation and
rotation. We propose an efficient and conservative particle-fluid coupling
scheme compatible with the FEM basis that enables convergence under grid
refinement of the two-way coupling terms. Efficient algorithms for neighbor
detection for particle-particle collision and particle-wall collisions are
adopted. The method is applied to a few different test cases and the results
are analyzed qualitatively. The results demonstrate the capabilities of the
implementation and the potential of the method for simulating large-scale
particle-laden flows in complex geometries.

##  1 Introduction

Fluid mechanics plays a crucial role in many physiological processes on health
and disease. Given recent advances in medical imaging, computational power,
and mathematical algorithms, real-time patient-specific computational fluid
dynamics is now becoming possible. We present a computational framework
capable of capturing fluid-particle interactions in complex medical-image
derived geometries that are much larger than the size of the particles and
with large number of particles. A four-way coupled volume-filtered Euler-
Lagrange solver with a stabilized finite element based incompressible flow
solver is developed.

Most of the development around particle-laden flows has been in the context of
industrial applications such as gas-solid flows in fluidized beds using
structured grids. In contrast, the application areas of interest here are
biological flows which often involve complex organic geometry. Finite element
method based incompressible flow solvers have been gaining popularity since
the development of stabilization schemes. They are attractive for their nice
mathematical properties and their ability to naturally deal with complex
geometries. FEM based solvers for Euler-Lagrange type particle-laden flows are
an active area of research. Casagrande et al. [4] have developed a stabilized
FEM based incompressible solver for particle-laden and porous flows that was
restricted to structured grids. Recently El Geitani et al. [7] have developed
a FEM based particle-laden flow solver for gas-particle flows in engineering
capable of handling unstructured hexahedral meshes for industrial
applications. Dufresne et al. [5] have implemented and demonstrated a
massively parallel particle-laden flow solver in the YALES2 code. This uses a
finite volume discretization for the fluid phase. While some “simple" complex
boundaries can be represented on structured grids through immersed boundary
methods for example, this is in general not very efficient and the fidelity of
the representation is limited by the grid size. Presence of large empty
regions in the domain can lead to excessive memory usage and computational
cost due to unused grid points. Unstructured grids are better because they
capture the complex geometry more accurately in comparison. But, unstructured
grids are more challenging to work with. The development of efficient
algorithms for particle tracking on unstructured grids is an active area of
research.

In this work we present the development of a versatile and massively parallel
framework for studying biological particles in subject-specific geometries by
combining (i) a recently-developed statistical hydrodynamic model for fluid-
particle flows; (ii) scalable Eulerian–Lagrangian algorithms; and (iii) state-
of-the-art parallel techniques for simulating biological flows in image-based
geometrical models. There are few salient aspects that set the current work
apart. The inter-phase coupling procedure developed here can be applied to a
wide range of particle and mesh sizes at once. It is shown to be highly
efficient, consistent and convergent. The particle-wall collision algorithm is
more efficient and avoid preprocessing steps. First we briefly describe the
mathematical model and particularly the volume-filtering procedure as applied
to the particle-fluid. Then we describe algorithmic developments that enable
scalable tracking of large number of particles within unstructured grids. Next
we describe the collision processing for inter-particle and particle-wall
collisions that can efficiently deal with complex non-convex boundary
collisions. Next, we describe a novel two stage particle-fluid interphase-
coupling transfer function that is efficient, conservative and convergent on
unstructured grids.

##  2 Volume-filtered Euler-Lagrange equations

The volume-filtered Euler-Lagrange approach to modeling particle-laden flows
is described in this section. The general idea is that the particle phase is
treated as being composed of Lagrangian particles with properties such has
diameter, density, velocity and potentially other scalar quantities like
temperature or species concentration. The effect of the fluid on the particles
is modeled using a so-called drag law, whereas the particle-particle
collisions are captured explicitly. This approach offers a good compromise
between expensive particle-resolved simulation methods and two-fluid/Eulerian-
Eulerian methods that rely heavily on sub-grid models to account for both
particle-particle and particle-fluid interactions. The equations for the fluid
phase and particle phase and the coupling is described next.

Figure 1: Schematic of Lagrangian particles inside a discretized fluid domain
along with a representation of the filtering kernel.

###  2.1 Fluid-phase equations

The equation for the fluid phase are obtained by applying a spatial filtering
operator on the incompressible Navier-Stokes equations, taking into account
the volume occupied by the particles and the momentum exchange between the
particles and fluid phase. The detailed derivation of the equations can be
found in [1, 3]. Figure 1 illustrates the filtering procedure on a system of
particles suspended in a fluid. The resulting continuity equation after
filtering reads

$\frac{\partial\phi_{f}}{\partial t}+\frac{\partial}{\partial
x_{i}}\left(\phi_{f}u_{i}\right)=0$ (1)

and the momentum equation is

$\rho\phi_{f}\frac{\partial u_{i}}{\partial t}+\rho\phi_{f}u_{j}\frac{\partial
u_{i}}{\partial x_{j}}=-\phi_{f}\frac{\partial p}{\partial
x_{i}}+\phi_{f}\frac{\partial\tau_{ij}}{\partial x_{j}}+\phi_{f}\rho
g_{i}+f^{p}_{i}$ (2)

where $\tau_{ij}=\mu\left(u_{i,j}+u_{j,i}\right)$, $\mu$ is the fluid
viscosity, $\rho_{f}$ is the fluid density, $g_{i}$ is the body force. The
term $f^{p}_{i}$ is the momentum feedback force from the particles acting on
the fluid. $\phi_{f}$ is the volume fraction of the fluid phase. The model
used to calculate the force on the particle due to the fluid is described
next.

###  2.2 Particle equations

The particles are modeled as rigid spheres. Their motion is calculated
according to Newtons second law for translation and rotation. For spherical
particles, due to the symmetry, tracking the absolute orientation is not
strictly necessary. The rotational velocity is required to calculate the
hydrodynamic torque based on the difference in the particle and fluid
rotational velocity. The translational motion for particle $p$ is given by

$m_{p}\odv{v_{i}}{t}=f_{i}^{h}+\sum_{a}f^{col}_{i,p\leftarrow
a}+(\rho_{p}-\rho_{f})V_{p}g_{i}$ (3)

where $m_{p}$ is the mass of the particle, $V_{p}$ is the volume of the
particle, $f_{i}^{h}$ is the hydrodynamic force exerted on the particle by the
fluid (elaborated further down). $f^{col}_{i,p\leftarrow a}$ is the force on
the particle due to collisions with other particles (and walls). Within the
summation, $p\leftarrow a$ denotes the collision of particle $p$ with particle
$a$. Gravity is included as a body force term denoted with the acceleration
due to gravity $g_{i}$. Since we do not include the resolved fluid stress
force in the particle momentum equation, we include a buoyancy force that
takes into account the difference in the density of the particle ($\rho_{p}$)
and the fluid ($\rho_{f}$). In the biofluids application area
$\rho_{p}/\rho_{f}\sim\mathcal{O}(1)$. If the resolved fluid stress force is
included in the particle momentum equation, then we would need to add the
total gravitational force on the particle, $\rho_{p}V_{p}g_{i}$ to the right
hand side of the particle momentum equation.

The angular momentum equation conservation for the _p_ th particle is given by

$I\odv{\omega_{i}}{t}=\tau^{d}_{i}+\sum_{a}\frac{d_{p}}{2}\epsilon_{ijk}n_{j}f^{tcol}_{k,p\leftarrow
a}$ (4)

where $I$ is the moment of inertia, $\omega_{i}$ is the angular velocity of
the particle. The first term on the right hand side is the hydrodynamic torque
due to drag and the second term the summation of the forces due to the
tangential component of the frictional force due to collisions. The torque
acting on a small isolated sphere rotating in a quiescent viscous fluid is
given by

$\tau^{d}_{i}=\pi\mu d^{3}_{p}\left(\omega^{f}_{i}-\omega^{p}_{i}\right)$ (5)

where $\omega^{f}_{i}$ is the angular velocity of the fluid interpolated to
the particle center and $\omega^{p}_{i}$ is the angular velocity of the
particle.

The hydrodynamic force can be written as a superposition of forces arising
from distinct mechanisms.

$f_{i}^{h}=f_{i}^{d}+f_{i}^{a}+f_{i}^{b}+f_{i}^{l}+f_{i}^{m}$ (6)

where $f_{i}^{d}$ is the quasi-steady drag force, $f_{i}^{a}$ is the added
mass force due to the acceleration of the fluid around the particle,
$f_{i}^{b}$ is the Basset history force to account for viscous effects due to
unsteady motion of the particle, $f_{i}^{l}$ is the Saffman lift force due to
the pressure distribution that develops on the surface of a particle in
velocity field with a non-zero gradient, and finally the Magnus lift force
$f^{m}$ due to the rotation of the particle. In general, this linear
superposition of individually identified hydrodynamic forces on the particle
is not well-founded, but is invariably used in literature. It can be shown to
hold in the low and high Reynolds number limits. In the absence of better
alternatives, we do the same and superpose the individually identified
hydrodynamic forces on the particle.

For the application areas of interest, Basset force and Magnus force are not
expected to be significant and are not included in the simulations. In some
Euler-Lagrange formulations in literature, there is an additional fluid stress
force (resolved stress term) added to the RHS. The particular Euler-Lagrange
formulation used here accounts for the resolved fluid force in the fluid
momentum equation and hence this term is not included in the particle momentum
equation. See Zwick and Balachandar [14] for a detailed discussion on the
different formulations and their equivalence.

The added mass force is given by

$f_{i}^{a}=\frac{1}{2}\rho V_{p}\left(\mdv{u_{i}}{t}-\odv{v_{i}}{t}\right)$
(7)

where $\mdv{~{}}{t}$ is the material derivative and $v_{i}$ is the velocity of
the particle. There are alternative forms in literature that account for the
finite size of the particle by adding a term called Faxen’s correction, that
has not been included here. $\odv{~{}}{t}$ is the time derivative in the
particle frame of reference. The Saffman lift force is calculated according to
the expression given by McLaughlin [8]. The equations are integrated using a
second order accurate scheme. The fluid solver is fully implicit and as a
result can take large time steps. If the fluid time step is larger than the
particle’s stable time step, the particle equations are sub-cycled in time to
avoid stability issues.

#### Quasi-steady drag force

The quasi-steady drag force is calculated using the relation provided by
Tavanashad et al. [11] for freely evolving suspensions of particles. The
correlation was developed from particle-resolved direct numerical simulations
conducted for a range of Reynolds numbers, volume fractions, and density
ratios. Most of the prior work on developing the drag relations assumed
particles to be much denser than the fluid phase, which is typically the case
in gas-solid flows relevant to many industrial applications such fluidized bed
reactors for example. Consequently, most of the the particle-resolved
simulations consider fixed particle packing and the particles are not allowed
to move. For the application being considered here, i.e., biological flows,
the difference in the density of the particles and the fluid is not very large
and the assumption of fixed particle packing might reduce the accuracy of the
drag relations. For $\rho_{p}/\rho_{f}<10$ the formula for the drag force is

$f^{d}_{i}\left(\phi,\operatorname{Re}_{m}\right)=f^{St}_{i}\left(1+0.15\operatorname{Re}_{m}^{0.687}\right)\left(78.96\phi_{p}^{3}-18.63\phi_{p}^{2}+9.845\phi_{p}+1\right)$
(8)

where $f^{St}_{i}=3\pi\mu d_{p}\left(u_{i}^{f}-v_{i}\right)$ is the Stokes
drag force and $\phi_{p}$ is the volume fraction of the particles. The
modified Reynolds number is defined as

$\operatorname{Re}_{m}=\frac{\rho_{f}d_{p}\phi_{f}\left|u_{i}^{f}-v_{i}\right|}{\mu}$
(9)

##  3 Numerics

###  3.1 Fluid discretization

The fluid phase equations equations are solved using the finite element
method. The filtered equations are written in their weak form and discretized
using the finite element method. The mass and momentum conservations can be
rearranged into a convenient form and written respectively as

$\displaystyle
u_{i,i}+\frac{1}{\phi_{f}}\left(\phi_{f,t}+u_{j}\phi_{f,j}\right)$
$\displaystyle=0,$ $\displaystyle\mathscr{L}_{i}=\rho u_{i,t}+\rho
u_{j}u_{i,j}+p_{,i}-\tau_{ij,j}-\frac{f^{p}_{i}}{\phi_{f}}$
$\displaystyle=0_{i}.$

The weak form of the above equations is obtained by multiplying by test
functions and integrating over the domain

$\displaystyle\int_{\Omega}\left\\{qu_{i,i}+\frac{q}{\phi_{f}}\left(\phi_{f,t}+u_{j}\phi_{f,j}\right)+w_{i}~{}\left(\rho
u_{i,t}+\rho
u_{j}u_{i,j}+p_{,i}-\tau_{ij,j}-\frac{1}{\phi_{f}}f^{p}_{i}\right)\right\\}~{}\mathrm{d}\Omega=0,$

where $q$ is the test function for the divergence constraint and $w_{i}$ is
the test function for the momentum equation. This is the standard (meaning
unstabilized) Galerkin weak formulation. We will refer to these terms as
Galerkin terms (to separate them from the stabilization terms that we will add
below). This is formulation is not stable in general. There are two issues to
overcome. One is the spurious upstream disturbance at high $Re$ and the other
is spurious modes in the pressure for certain combination of velocity-pressure
elements related to the LBB or inf-sup condition. Same type of elements for
velocity and pressure in particular do not work. There is some attraction to
use same element types for both velocity and pressure. To remedy the high $Re$
problem and to remove restrictions on the choice solution function spaces we
add stabilization terms to the formulation. Let us denote the last LHS by
$B_{G}(w_{i},q;u_{i},p)$. With the stabilization terms added the formulation
is

$\displaystyle B(w_{i},q;~{}u_{i},p)$ $\displaystyle=0$ $\displaystyle
B(w_{i},q;~{}u_{i},p)$ $\displaystyle=B_{G}(w_{i},q;~{}u_{i},p)$
$\displaystyle+\sum_{e=1}^{n_{el}}\int_{\bar{\Omega}_{e}}\left\\{\tau_{M}(u_{j}w_{i,j}+\frac{q_{i}}{\rho})\mathscr{L}_{i}+\tau_{C}w_{i,i}u_{j,j}\right\\}~{}\mathrm{d}\bar{\Omega}_{e}$
$\displaystyle+\sum_{e=1}^{n_{el}}\int_{\bar{\Omega}_{e}}\left\\{w_{i}\rho\overset{\Delta}{u_{j}}u_{i,j}+\bar{\tau}\left\\{\frac{\rho\overset{\Delta}{u_{j}}}{\tau_{M}}\right\\}w_{i,j}\left\\{\frac{\rho\overset{\Delta}{u_{k}}}{\tau_{M}}\right\\}u_{i,k}\right\\}~{}\mathrm{d}\bar{\Omega}_{e},$

The details of the stabilization scheme and the complete details of the
numerical approach can be found in [12]. In brief, consistent terms are added
to the weak form to stabilize the system. In the continuous limit, the
stabilization terms vanish and we are left with a consistent set of equations.
The non-linear fluid equations are solved implicitly in a monolithic fashion
using a Newton-Raphson method. The time integration scheme used is second
order generalized-$\alpha$ method suitably modified for the system of
differential algebraic equations that arise from the discretization of the
fluid equations. We use linear tetrahedral elements for all variables.

The implementation is done within the CRIMSON framework [2]. CRIMSON
specializes in simulation of blood flow in the vasculature. It has advanced
capabilities to read-in medical images and allows the user to create a
corresponding geometry and mesh and apply complex physiologically accurate
dynamic boundary conditions. The underlying flowsolver is a general purpose
finite element solver, PHASTA [12] that solves the incompressible Navier-
Stokes equations on unstructured grids. The solver is parallelized using MPI
and has been demonstrated to scale to 100,000 cores[13].

###  3.2 Particle discretization

The particle related computations are described in this section. First we
describe how the particles are tracked through the domain as they are carried
by the flow. Then we describe how the collisions are detected and processed.

####  3.2.1 Particle tracking

Tracking particles in an unstructured grid is a challenging task. The main
difficulty is that the elements are not arranged in a regular fashion. So, we
cannot use a simple indexing scheme to identify the element that contains the
particle. During the course of the simulation, as teh particles are advected
by the flow, the particles will move from one element to another. After each
particle has been advected, we identify the new element host for each particle
by executing a delaunay search with the last known element as the starting
point. For the applications of interest, the time-step is small enough that
the particles do not move very far and the delaunay search works well.

When the particles are initialized, as a preprocessing step we identify and
record the starting element for each particle to enable the quick delaunay
search.

The particle tracking and collision detection is implemented using a two stage
approach. The first stage is to identify the element that contains the
particle. The second stage is to localize the particle within the element. The
first stage is implemented using a cell list method. The second stage is
implemented using a delaunay search method. In prior work Rydquist and Esmaily
[10] have used a cell list method for particle tracking. The idea is to
overlay a Cartesian grid on top of the unstructured mesh and assign each
cartesian cell with all the unstructured mesh elements that intersect with
this cell.

First, we initialize a integer array for the background grid. We will store
just one value per cell. This can be a hashmap if the domain has a lot of
empty space and memory usage is a concern. Lets say that the side of the cell
is approximately twice as large as a characteristic length of the largest
element in the unstructured mesh.

Then, we loop over the elements and for each element find indices of the cell
that contains the first vertex of this element. If this cell has not already
been assigned, assign this element id to this cell. At this point one might
think that that all relevant cells have been assigned a value. If we draw the
cells that have been assigned overlaid on the mesh, we would like to see that
the cartesian cells cover the entire region. But it is very likely that there
are cells which did not have the vertex that was tested( or any vertex) inside
them. So, there will be "holes" if the cartesian cells are too small.
Additionally, near the boundaries (particularly curved or inclined) it is
possible that there are parts of the element that do not have a cartesian cell
covering it. If the cartesian cells are not too small ( i.e. large enough that
each cell contains at least one whole element), then the holes do not appear.
To deal with the potential lack of coverage with boundary elements, the
following is done. We loop over the cartesian cells and for each cell check if
the immediate neighbors are assigned. if not we make a note of the indices and
the value that will be assigned to this adjacent cell will be the value
assigned this cartesian cell. This needs to be done in two stages,
identification and assignment, because we wouldn’t know if a cell was assigned
the first time or during this step. The effect will be that we are padding the
cartesian cells with one layer of cells. This ensures full coverage of
boundary elements.

####  3.2.2 Collision processing

This section described the algorithms adopted for scalable calculation of
particle-particle and particle-wall collisions are described. It is well known
that a naive implementation of nearest neighbor search for particle-particle
collisions scales as $\mathcal{O}(N_{p}^{2})$. Similarly a naive
implementation of particle-wall collisions scales as $\mathcal{O}(N_{p}N_{w})$
where $N_{w}$ is the number of wall elements. This is a non-starter even for a
moderate number of particles. To enable simulation of $\mathcal{O}(10^{9})$
particles we need to avoid the quadratic and multiplicative scaling of the
collision detection algorithms. Such algorithms for both particle-particle and
particle-wall collisions are described below.

##### 3.2.2.1 Particle-particle collisions

Figure 2: Using a distance field to the partition boundary, we identify
particles on neighboring processes that might collide with particles on a
given process and communicate these ghost particles, the particles that are in
the green shaded region, to the process that owns the blue shaded partition.

The flow solver uses domain decomposition to solve the equations in parallel.
Working within this framework, we store the particles on the process that
contains the element that contains the center of the particle. A simple
example of a meshed domain that has been decomposed into three partitions is
shown in fig. 2. Let us consider the partition highlighted in blue. To capture
all the collisions of the particles in this partition, the process that owns
this partition needs to know about particles on neighboring processes that are
nearby the shared partition boundary within some cut-off distance. This
distance will depend on the size of the particles. As we use a soft sphere
model for collisions, the cut-off distance should be at least the diameter of
the largest particle. If there are some guarantees about the relative size of
the particles with respect to the mesh size, then, the elements in the cut off
region could be counted based on their graph distance, i.e. we could determine
this halo zone in terms of layers of cells. But in general such assumptions
are not possible and a general approach is called for. Here, we identify the
the elements in the cut off region by calculating the distance field from the
from the partition boundary to the interior of the partition. Then, the
particles in the cut off region are identified by interpolating the distance
field to the particle centers. These are then communicated to the
corresponding neighboring process. This is illustrated in fig. 2 where the
elements highlighted in green are the elements that contain particles that can
potentially collide with particles in the blue partition.

Figure 3: After ghost particles are communicated, a cell list is used to
identify nearby particles for collision detection. To identify potential
collisions for particles in the green cell, we only need to check the
particles in the green cell and the neighboring blue-outlined cells.

Once all the particles are available on a process, the remaining work is
processor local and does not depend on parallelization details. Using a cell
list we divide the bounding box containing the particles into Cartesian cubes.
The size of the cube is based on the largest particle diameter. We loop over
the particles and assign them to the cell that contains the particle center.
Once this is done for all particles, to identify collisions for a given
particle, we only need to check the particles in the cell that contains the
particle and the neighboring cells ($1+26=27$ neighbors in 3D). This is
illustrated in fig. 3 in two dimensions. To compute all possible collisions
for a particle in the green cell, we only need to check the particles in the
green cell and the neighboring blue-outlined cells. This is a significant
reduction in the number of particles that need to be checked for collisions.
The scaling of the algorithm with the number of particles is shown in fig. 4.
The choice of the size of the cell list box is important. For optimal results,
Rydquist and Esmaily [10] has shown that the number of cells should
approximately be equal to the number of particles.

Figure 4: Scaling of the collision detection algorithm with number of
particles

Since we are dealing with complex tortuous geometries, the distribution of the
particles within the bounding box will not be uniform. This can lead to
excessive memory usage if we allocate a 3d array of lists to store the
particles in each cell. To avoid this issue, we use a locality sensitive
hashing scheme to store the particles in a unsorted key-value store. The key
is a hash of the cell indices and the value is a list of particles. This is
implemented using the std::unordered_map data structure in C++, also known as
a hash map. With this method, if a cell is empty, we do not need to allocate
any memory for it. This is a significant advantage over the 3d array of lists
approach. The particular hash function used is the 3D Morton index or Z-order
curve. This is a space filling curve that reversibly maps a 3D point to a 1D
index. The mapping is such that points that are close in 3D space are mapped
to indices that are close in 1D space. This is a desirable property for the
hash function because it means that the particles that are close in 3D space
will be mapped to keys that are close in 1D space. This will lead to a more
uniform distribution of the particles in the key-value store. In addition,
since this is a reversible hash, by knowing the key, we can recover the 3D
cell index for any particle. The scaling of this collision detection algorithm
with the number of particles is shown in fig. 4 going up to a million
particles. We get the optimal linear scaling with the number of particles.

##### 3.2.2.2 Particle-wall collisions

The collision detection can be accelerated if we only check the boundary
collisions of only those particles that are near the boundaries and ignore the
particles in the interior. To enable quick identification of near wall
particles, we initialize a distance field due to the boundary of the domain in
the interior of the domain. At any point in the interior, the distance field
is the shortest distance to the nearest boundary. This is illustrated in fig.
5 for a 2D aorta geometry. Near wall particles are identified by interpolating
the distance field to the particle center and checking if the distance is less
than the about 1.2 times radius (including a factor of safety). Once the near
boundary particles have been flagged, we use a another cell-list like
structure and assign the boundary triangles to the cells. This reduces the
number of boundary features that need to be checked to only those that are
nearby the particle and skips triangles, edges and vertices that are far away.
Then, for each particle we easily identify the cell and test for collision
with the boundary faces, edges and vertices in the neighboring cells and
calculate the collision force. We first calculate a distance field in the
interior of the domain. To quickly eliminate particles in the interior, we use
the FEM interpolating functions to calculate the distance at the particle
centers.

Figure 5: Contour plot of the boundary distance field for a 2D aorta geometry.
The farther away from the boundary the darker the color.

In general, the particles near the domain boundary could collide with a
boundary feature such as a face, an edge, a vertex or any permutation of
multiple of them. Robust and unambiguous processing of collisions is made
possible by using the concept of Voronoi regions[5]. The near boundary region
is divided into regions based on which boundary feature is nearest to the
point in the region. This is illustrated in fig. 6. Particles near the
boundary are assigned to the region corresponding to the feature that is
nearest to the particle center. Collision force is calculated only if the
particle is in the Voronoi region of the boundary feature it is intersecting
with. The algorithm adopted here first checks for face collisions. These are
always valid and are registered immediately. For any faces that are
registered, we void the edges and vertices. That is, we ignore any
intersection of this particle with any edge or vertex that were part of
previously registered faces. Then, we do a pass over the edges and register
edge collisions and void the vertices for the edges that are registered.
Finally, we do a pass over the vertices and register vertex collisions.

Figure 6: Voronoi region near the domain boundary for three different types of
geometric features that make up the boundary, i.e., faces, edges and vertices.

The approach presented here deals with correctly with all possible collision
scenarios regardless of the local curvature characteristics of the boundary
surface. It is a mathematically exact and exhaustive procedure without any
user-specified tolerance that might require tuning. In literature, the idea of
Voronoi region has been used but it required an expensive preprocessing step
to actually build and store the Voronoi information. The algorithm described
here implicitly implements the Voronoi region idea and does not require a
preprocessing step and has improved efficiency while retaining the robustness
and unambiguous nature of the original Voronoi region approach.

To calculate the actual collision force, once the collision normal and
penetration depth are calculated, we treat the wall as a particle with
infinite mass and calculate the collision force as if it were a particle-
particle collision using the same soft sphere model.

##  4 Particle-fluid coupling

The coupling between particles and fluid is a crucial aspect of Euler-Lagrange
simulations. Depending on the specific characteristics of the particle-laden
flow being simulated, different levels of physical interactions between the
particle and fluid phases can be incorporated in the simulation. In one-way
coupling, the fluid exerts a drag force on the particles, but the particles do
not affect the fluid. This is typically observed in dilute flows where the
particle volume fraction is small and the particles have minimal impact on the
fluid flow. Two-way coupling involves interactions in both directions, where
the particles affect the fluid and the fluid affects the particles. This is
commonly seen in dense flows where the particles significantly influence the
fluid flow and vice versa. Four-way coupling takes into account not only the
interactions between particles and the fluid, but also considers inter-
particle collisions. This is suitable for dense flows with strong particle-
particle and particle-wall collisions. Within two-way and four-way coupling,
it is possible to ignore the excluded volume effects of the particles if the
particle volume fraction is small, but still consider the momentum exchange
between the particles and the fluid.

In one-way coupled simulations, we need to compute the fluid velocity,
vorticity, fluid stresses etc. at the particle center to be able to calculate
all the hydrodynamic forces on the particle and integrate the particle
equations. For higher levels of coupling in addition to the Eulerian to
Lagrangian transfer, we also need to have Eulerian representations of particle
centered quantities, mainly the volume fraction field and feedback force.

Ideally, transfers in both directions would be done consistently using the
same filtering kernel function by identifying the support of the kernel
function and integrating over the support. But, this is not practical due to
the $\mathcal{O}(N_{p}N_{sn})$ complexity involved where $N_{sn}$ is the
number of nodes in the support of the kernel function. The time complexity
further hides the fact that identifying the support nodes in a unstructured
grid is a non-trivial task in itself. In literature, there have been a variety
of ad hoc methods that have been used to project particle centered quantities
to the Eulerian mesh. But these are lacking in either accuracy, consistency or
are grid-dependent. These properties become important when dealing with
particle sizes about equal to or greater than the local mesh size and can
severely compromise the accuracy of the simulation. In order to sidestep the
complexity issues, various alternative methods of transferring information
between the particle and fluid phases have been proposed in literature. There
are different considerations for the Eulerian-to-Lagrangian and Lagrangian-to-
Eulerian transfers. There are no conservation and boundedness issues when
interpolating from the Eulerian mesh to the particles. It is a straight-
forward interpolation and one can devise efficient schemes to do this. The
only consideration is the accuracy of the interpolation.

###  4.1 Particle-to-fluid transfer

The Lagrangian-to-Eulerian transfer needs more care. For example, a particle-
in-cell type approach would assign the whole volume of the particle to the
Eulerian cell it is contained in. If we project using a mesh-size dependent
kernel, we will end up with a mesh dependent result that can give rise to
spurious and unbounded values. The volume fraction field in particular needs
to lie in $[0,1]$ to be physically meaningful. The random close packing limit
for mono-disperse spheres is $\sim 0.64$. So, in reality in a non-jammed
particle-laden flow the particle volume fraction will be in $[0,0.64]$.

In literature, there have been a variety of methods that have been used to
project particle centered quantities to the Eulerian mesh. But these are
lacking in either accuracy, consistency or are grid-dependent. For instance,
Casagrande et al. [4] use a particle-in-cell type approach where the volume
occupied by a particle is assigned to the element containing the center. There
are no conservation issues and this works fine if the particles are much
smaller than the mesh size. If not, there is a possibility of unphysical
values appearing in the projected quantities. In [7] a similar idea is used to
get element-constant values that are then projected to nodal values while
artificially restricting the nodal values to be within physical limits.
Obtaining a volume fraction field from a particle configuration was
exclusively considered in [6], where different methods are reviewed.

Figure 7: Lagrangian to Eulerian projection of the volume fraction for one
particle centered at the origin in a 1D mesh using the FEM interpolating
polynomials. As the mesh is refined, the projected value diverges

Pusztay et al. [9] use the finite element basis functions to perform the
Lagrangian to Eulerian projection and directly obtain nodal values on the
mesh. The process can be formulated as enforcing a weak equivalence between
the particle volumes as a collection of $\delta$-functions at the particle
centers and the Eulerian mesh representation of the particle volume fraction
field. This is conservative, but it has same boundedness and convergence
issues as the PIC-type method. In addition, depending on the specifics of the
quadrature rules used, this method can generate oscillations and potentially
negative values locally. The divergent result of this method is demonstrated
in fig. 7. It shows the projected volume fraction field on a 1D finite element
mesh due to one particle at the origin under mesh refinement. As the mesh is
refined, the projected value diverges.

###  4.2 Two stage projection

Here we outline a novel two-stage method that does not suffer from the issues
described before and is efficient and scalable. It is a generalization of the
two stage method described in Capecelatro and Desjardins [3] to unstructured,
anisotropic and inhomogeneous meshes. The first step is along the lines of
[9], with one modification to ensure that we preserve the original sign of the
quantity being projected. Let the particle distribution be denoted by

$\psi=\sum\limits_{i}^{N_{p}}V_{i}\delta(\bm{x}-\bm{x}_{i})$

where the $V_{i}$ are the particle volumes and $\bm{x}_{i}$ are the particle
centers. The finite element representation of the volume fraction field cane
be written as

$\phi_{p}=\sum\limits_{i}^{N_{n}}\varphi_{i}N_{i}(\bm{x})$

where $N_{i}$ are the finite element basis functions and $\varphi_{i}$ are the
nodal values. The first step is to enforce the equivalence between the
particle distribution and the finite element representation of the volume
fraction field in a weak sense. This is done by enforcing the following

$\int_{\Omega}\phi_{p}N_{i}(\bm{x})~{}\mathrm{d}\Omega=\int_{\Omega}\psi
N_{i}(\bm{x})\mathrm{d}\Omega$ (10)

This can be simplified and written as a linear system of equations of the form

$[M]~{}\\{\varphi\\}=[V]~{}\\{w\\}$ (11)

where $[M]$ is the mass matrix, $\\{\varphi\\}$ is the vector of nodal values
and $\\{w\\}$ is the vector of particle volumes. $[V]$ is a sparse matrix
whose rows correspond to the particle centers and columns correspond to the
mesh nodes. The mass matrix in eq. 11 should be lumped for a couple of
different reasons. First, lumping avoids oscillatory and negative values in
the projected field. It also diagonalizes the linear system and the solution
can be obtained by simple division at each node. This result is still mesh
dependent and cannot in general be used if the particle sizes are comparable
or larger than the mesh size.

Figure 8: Same case as fig. 7 projected using the two-stage method. Under mesh
refinement, the result converges to a Gaussian of the prescribed width (i.e.
filter size $\delta_{f}$). Right shows rate of convergence

The second step is to smooth this projected field using a diffusion equation,
using a carefully chosen diffusivity that takes into account the local mesh
size to give a mesh independent smoothing effect. The diffusion equation is
discretized using the same finite elements and advanced using pseudo-time-
stepping on the same domain as the fluid. This can be written as

$\pdv{\phi_{p}}{\tau}=\nabla\cdot\left(\mathscr{D}\nabla\phi_{p}\right)$ (12)

with a diffusivity defined as

$\mathscr{D}=\max{\left(\frac{\delta_{f}^{2}-\Delta x^{2}}{16\ln 2},0\right)}$
(13)

where $\delta_{f}$ is the filter size/smoothing-length-scale set to be
$4d_{p}$ and $\Delta x$ is a nominal local mesh size that is calculated as the
average of the mesh edge lengths for each node. Here, $\tau$ is the pseudo
time variable and eq. 12 is advanced from $\tau=0$ to $\tau=1$. This
corresponds to a smoothing length scale of $\sqrt{\mathscr{D}\tau}$ which
corresponds to a Gaussian with a full width at half maximum of $\delta_{f}$
i.e. it is equivalent to analytically transforming a Dirac delta function to a
Gaussian of width $\delta_{f}$ while conserving the integrated value. While
the method is demonstrated here on a 1D case, the extension to higher
dimensions is quite straight forward and we have verified that the 3D version
retains the desired properties.

##  5 Example cases

The method developed here is applied to a few different particle-laden flow
cases. This is preliminary work and the focus is on demonstrating the
capabilities of the implementation. The results presented here are presently
only qualitative. In addition to a simple pipe geometry, we demonstrate the
capabilities of the method by simulating particles flowing through a
simplified bifurcation geometry that one might find in the arterial
vasculature.

###  5.1 Pressure driven flow in a pipe

The first case considered is pressure driven flow in a cylindrical tube. The
geometry is shown in figure 9. The Reynolds number based on the pipe diameter
is 100. There are 10,000 particles are initialized in as a spherical bolus
near the inlet boundary. The domain is discretized using 24,000 elements. The
inlet velocity is specified as parabolic flow profile. The outlet is a zero
pressure boundary condition and the walls are no-slip boundaries. The particle
are denser than the carrier fluid, $\rho_{p}=5\rho_{f}$ and gravity is
included. This case is solved in parallel with 4 MPI processes.

Figure 9: Fully coupled simulation of particles in a pressure driven pipe
flow. Results at different instances are shown. In each frame, the figure on
the left shows contour plot of the volume fraction of the particles on a
vertical slice through the center of the cylinder. The figure on the right
shows a contour plot of the magnitude of the velocity field along with the
particles themselves.

In figure 9 we see that the particles are transported down the cylinder by the
fluid. Due to the presence of gravity particles also settle along the bottom.
The fully-coupled nature of the simulation is evident from the velocity
contours. Initially, the particles are at rest and block the fluid. As a
result, he fluid has to accelerate around the particles. This delays the
development of the parabolic profile downstream of the particles. The
particles slowly accelerate due to the fluid and gradually pick up speed and
are carried out of the cylinder through the outlet.

###  5.2 Flow in bifurcation geometry

#### Flow in a simple bifurcation

Another example demonstrated here is particle-laden flow through a bifurcation
geometry shown in fig. 10. The fluid domain is vessel that splits into two
smaller vessels, that resembles a typical feature of the arterial vasculature.
The vessel geometry is discretized using 220K tetrahedral elements. We
initialize 100,000 particles in the inlet in a cylindrical bolus with a mean
particle volume fraction in the bolus of $\sim 0.35$. The particles and fluid
are initially at rest. The inlet velocity is specified with a parabolic
profile corresponding to a Reynolds number of 250 based on the inlet flow and
diameter. There is no gravity and the particles density is twice that of the
fluid. The outlet vessels are coupled to a 0D lumped parameter model that
models the downstream vasculature. The lumped parameter model is a 3-element
Windkessel model. Since this is a steady flow simulation, after the lumped
parameter capacitor is charged, the outflow boundary conditions reduce to
outlets with a resistance.

Figure 10 shows the initial condition in the top left panel and snapshots at
two subsequent instances as the simulation progresses. Within each subfigure,
the left panel shows a contour plot of the velocity field magnitude along with
the particle positions. The right panel is a close-up view of the bifurcation
region looking up stream into the inlet vessel. As the particles come to the
bifurcation, they collide with the boundary walls and get split into both the
downstream vessels following complex inertial trajectories.

Figure 10: Fully coupled simulation of particles flowing through a
bifurcation. The initial condition is shown in (a). Snapshots of the velocity
contour and particle locations at subsequent times are are shown in (b) and
(c)

#### Flow in a bifurcation with a stenosis

The final example shown is a modification of the bifurcation geometry from the
previous example. Here, we introduce a stenosis in one of the daughter
vessels, to simulate a narrowing of the vessel due to atherosclerosis. The
rest of the conditions are unchanged. The results for this case are shown in
fig. 11.

Figure 11:  Comparison of particles flowing through a bifurcation with and
without a stenosis. In (d) the particles can be seen to be flowing back away
from the outlet in the stenosed branch. This is due to the formation of a
circulating flow topology caused by the stenosis.

Within each subfigure in fig. 11, the previous healthy geometry is on the
right and the stenosed geometry is on the left. We compare these two
simulations at same instances in time. Till subfigure (b) the the particle
configurations do not show any significant difference between the two cases.
In fig. 11c we see that the particles are accelerated through the narrowing
and reach the outlet much quicker than before. In the final subfigure, in the
stenosed geometry we see a particles form a spur-like shape and flowing
backward. This is due to a recirculation zone that forms downstream of the
stenosis, a well known feature of flow through stenosed vessels.

These examples are meant to be a qualitative demonstration of the capabilities
of the method. In the future, this work will be followed up with validation
and performance evaluation and scaling studies. We will also apply this method
to more larger and complex geometries derived from patient-specific medical
images.

##  6 Conclusions

To summarize, we have presented a computational framework for modeling large
scale particle-laden flows in complex geometries that enable subject-specific
analysis. The framework is based on a volume-filtered Euler-Lagrange approach
that uses a finite element method for the fluid phase and a discrete element
method for the particle phase. The fluid phase is solved on a three-
dimensional unstructured grid using a stabilized finite element method. The
particle phase is modeled as rigid spheres and their motion is calculated
according to Newton’s second law for translation and rotation. The
hydrodynamic force on the particles is calculated using a recently developed
correlation for freely evolving suspensions of particles. The method is
applied to a few different particle-laden flow cases. The results are of a
qualitative nature and are not meant to be quantitative. The results
demonstrate the capabilities of the implementation and the potential of the
method for simulating large-scale particle-laden flows in complex geometries.
We intend to follow up this work with validation and performance evaluation
and scaling studies. We will also apply this method to more larger and complex
geometries derived from patient-specific medical images and other biofluids
applications.

## References

  * Anderson and Jackson [1967] T. B. Anderson and Roy Jackson.  Fluid Mechanical Description of Fluidized Beds. Equations of Motion.  _Industrial & Engineering Chemistry Fundamentals_, 6(4):527–539, November 1967.
  * Arthurs et al. [2021] Christopher J. Arthurs, Rostislav Khlebnikov, Alex Melville, Marija Marčan, Alberto Gomez, Desmond Dillon-Murphy, Federica Cuomo, Miguel Silva Vieira, Jonas Schollenberger, Sabrina R. Lynch, Christopher Tossas-Betancourt, Kritika Iyer, Sara Hopper, Elizabeth Livingston, Pouya Youssefi, Alia Noorani, Sabrina Ben Ahmed, Foeke J. H. Nauta, Theodorus M. J. Van Bakel, Yunus Ahmed, Petrus A. J. Van Bakel, Jonathan Mynard, Paolo Di Achille, Hamid Gharahi, Kevin D. Lau, Vasilina Filonova, Miquel Aguirre, Nitesh Nama, Nan Xiao, Seungik Baek, Krishna Garikipati, Onkar Sahni, David Nordsletten, and C. Alberto Figueroa.  Crimson: An open-source software framework for cardiovascular integrated modelling and simulation.  _PLOS Computational Biology_ , 17(5):e1008881, May 2021.
  * Capecelatro and Desjardins [2013] Jesse Capecelatro and Olivier Desjardins.  An Euler-Lagrange strategy for simulating particle-laden flows.  _Journal of Computational Physics_ , 238:1–31, 2013.
  * Casagrande et al. [2017] Marcus V. S. Casagrande, José L. D. Alves, Carlos E. Silva, Fábio T. Alves, Renato N. Elias, and Alvaro L. G. A. Coutinho.  A hybrid FEM-DEM approach to the simulation of fluid flow laden with many particles.  _Computational Particle Mechanics_ , 4(2):213–227, April 2017.
  * Dufresne et al. [2020] Yann Dufresne, Vincent Moureau, Ghislain Lartigue, and Olivier Simonin.  A massively parallel CFD/DEM approach for reactive gas-solid flows in complex geometries using unstructured meshes.  _Computers & Fluids_, 198:104402, February 2020.
  * El Geitani and Blais [2023] Toni El Geitani and Bruno Blais.  Quadrature-Centered Averaging Scheme for Accurate and Continuous Void Fraction Calculation in Computational Fluid Dynamics–Discrete Element Method Simulations.  _Industrial & Engineering Chemistry Research_, 62(12):5394–5407, March 2023.
  * El Geitani et al. [2023] Toni El Geitani, Shahab Golshan, and Bruno Blais.  Toward High-Order CFD-DEM: Development and Validation.  _Industrial & Engineering Chemistry Research_, 62(2):1141–1159, January 2023.
  * McLaughlin [1993] John B. McLaughlin.  The lift on a small sphere in wall-bounded linear shear flows.  _Journal of Fluid Mechanics_ , 246:249–265, January 1993\.
  * Pusztay et al. [2022] Joseph V. Pusztay, Matthew G. Knepley, and Mark F. Adams.  Conservative Projection Between Finite Element and Particle Bases.  _SIAM Journal on Scientific Computing_ , 44(4):C310–C319, August 2022.
  * Rydquist and Esmaily [2020] Grant Rydquist and Mahdi Esmaily.  An optimal O(N) scheme for simulations of colliding, particle-laden flows on unstructured grids.  _Journal of Computational Physics_ , 420, November 2020.
  * Tavanashad et al. [2021] Vahid Tavanashad, Alberto Passalacqua, and Shankar Subramaniam.  Particle-resolved simulation of freely evolving particle suspensions: Flow physics and modeling.  _International Journal of Multiphase Flow_ , 135:103533, February 2021.
  * Whiting and Jansen [2001] Christian H. Whiting and Kenneth E. Jansen.  A stabilized finite element method for the incompressible Navier-Stokes equations using a hierarchical basis.  _International Journal for Numerical Methods in Fluids_ , 35(1):93–116, January 2001.
  * Zhou et al. [2010] Min Zhou, Onkar Sahni, H. Jin Kim, C. Alberto Figueroa, Charles A. Taylor, Mark S. Shephard, and Kenneth E. Jansen.  Cardiovascular flow simulation at extreme scale.  _Computational Mechanics_ , 46(1):71–82, June 2010.
  * Zwick and Balachandar [2020] David Zwick and S. Balachandar.  A scalable Euler–Lagrange approach for multiphase flow simulation on spectral elements.  _International Journal of High Performance Computing Applications_ , 34:316–339, May 2020.