11institutetext: Delft University of Technology, Delft, The Netherlands
11email<EMAIL_ADDRESS>
11email<EMAIL_ADDRESS>
# Parameterizing Federated Continual Learning for Reproducible Research

Bart Cox  0000-0001-5209-6161 Jeroen Galjaard  0000-0003-3681-7226 Aditya
Shankar  0009-0009-3046-8724 Jérémie Decouchant  0000-0001-9143-3984 Lydia Y.
Chen  0000-0002-4228-6735

###### Abstract

Federated Learning (FL) systems evolve in heterogeneous and ever-evolving
environments that challenge their performance. Under real deployments, the
learning tasks of clients can also evolve with time, which calls for the
integration of methodologies such as Continual Learning. To enable research
reproducibility, we propose a set of experimental best practices that
precisely capture and emulate complex learning scenarios. Our framework,
Freddie, is the first entirely configurable framework for Federated Continual
Learning (FCL), and it can be seamlessly deployed on a large number of
machines thanks to the use of Kubernetes and containerization. We demonstrate
the effectiveness of Freddie on two use cases, (i) large-scale FL on CIFAR100
and (ii) heterogeneous task sequence on FCL, which highlight unaddressed
performance challenges in FCL scenarios.

###### Keywords:

Federated Continual Learning Resource and Data heterogeneity

##  1 Introduction

Federated Learning (FL) [12] performs distributed optimization thanks to a
central federator server that maintains a global model using model updates
computed by clients. It is common for data to be distributed among the clients
of an FL system in a non-independent and identically distributed (non-IID)
way. Moreover, in practice, client learning tasks also evolve over time.
Continual Learning (CL) [5] is a technique that addresses the scenario where a
model is continuously trained on evolving client tasks.

One of the key challenges in CL is catastrophic forgetting: parameters or
semantic representations learned for past tasks drift under the influence of
new tasks. Three categories of techniques address this challenge [6]. Replay
mechanisms, like GEM [8] and DGR [15], retain or generate data from earlier
tasks for new task adaptation, which allow the network to revise previously
learned tasks. Regularization techniques, such as EWC [5], penalize the
divergence of model parameters, preventing the adaptation process on new tasks
from deviating too far from the model learned on prior tasks. Parameter
isolation methods use specific weights of the network for the task at hand,
i.e., use a mask to freeze the weights of other tasks [10].

Continual Learning allows a client to learn from its previous tasks if
features are repeated over time. Federated Continual Learning [18] (FCL)
combines CL and FL, enabling clients to indirectly learn from each other.
Existing CL frameworks do not take this indirect learning into account and
therefore provide limited support for Federated Continual Learning.

In addition, reproducing FCL results that were obtained in deployment is
difficult. For example, experimental environments are often tightly controlled
and steady, while real-world environments are often dynamic and heterogeneous.
In addition, clients might be punctually busy processing co-located tasks.
Several FL simulation [9, 14] and emulation [1, 14] frameworks have been
proposed, but they cannot be easily extended to support heterogeneous data,
learning tasks and hardware platforms. In addition, frameworks that focus on
enabling large-scale FL experiments impose a significant overhead to manage
the execution or require the use of a strict pipeline.

In this paper, we address the lack of a scalable yet flexible framework for
reproducible FCL experiments. Overall, we make the following contributions:

  * •

We identify key requirements for FL and FCL emulation: ease of use,
reproducibility, support for complex workloads, and resource heterogeneity.

  * •

We develop Freddie—a framework for Federated and distributed machine-
learning—the first open
source111https://gitlab.ewi.tudelft.nl/dmls/publications/freddie framework
that addresses these requirements. Freddie supports small scale deployments,
i.e., single machine simulations, and large scale emulation over self-managed
and cloud systems using Kubernetes. Freddie enables the emulation of both data
and resource heterogeneity.

  * •

We provide benchmarking generating methods for FCL that explore both data and
task heterogeneity across clients —realistic workloads tailored for Federated
Continual Learning systems.

##  2 Related Work

Federated Learning (FL). Existing FL frameworks support a fixed set of
learning tasks across clients. Flower [1] provides a client-server framework
that needs to be manually started on different devices. Differently, Fate [3]
focuses on providing a secure and production-ready federated learning setup.
Fate supports Kubernetes deployments but requires the use of its pipelines to
run experiments. Although this might provide desirable security additions for
production, it also tends less to prototyping and active research needs.
Besides research endeavours, popular deep learning platforms such as TorchX
and Tensorflow Federated can respectively run distributed and federated
experiments at scale, but they lack the flexibility to use other ML libraries.

Continual Learning (CL). FACIL [11], PyCIL [19], and Pycontinual [4] provide
CL frameworks and CL algorithms such as LwF, iCaRL, EWC, and GEM. Continual
World [17] adds a simulation world for robotics tasks for Continual
Reinforcement Learning. Avalanche [7] is focused on reproducible End-to-End
Continual Learning. The aforementioned frameworks support CL only on a single
machine. FedWEIT [18] combines parameter isolation and regularization and
extends CL to a federated setting. However, it does not consider the impact of
task sequences on the global model’s quality. Last but not least, current FL
frameworks cannot be easily extended to support CL scenarios where the output
types evolve.

Continual Learning methods leverage task IDs during training and evaluation,
enabling the exploitation of a specific set of model weights or restricting
the output classes based on a task’s ID. This is also known as Task-
Interactive Learning (_task-IL_) and Domain-Interactive Learning (_domain-IL_)
[16], which are both supported by Freddie.

##  3 Specifying Data Heterogeneity and Learning Workloads

This section briefly surveys system requirements addressed by Freddie.

Classical FL parameters. The number of clients, the federators’ aggregation
and client selection strategies must be configurable. In addition, FL-related
and common hyperparameters, such as the training epochs, learning rates, etc.,
must be configurable.

Statistical heterogeneity. As in FL, local data distributions remain of high
importance, and their non-iidness should be configurable. In the context of
FCL, local distributions also limit the tasks that clients might be able to
train for.

Resource heterogeneity. It should be possible to specify the computing power
of the clients and of the federator, and the characteristics (latency,
throughput) of the network links that interconnect them.

Task description. The task sequence of each client can be specified. Non-IID
task distribution can be assimilated to the situation where clients learn
tasks with high intra-task variance, e.g. due to different domains. In such
settings, it is often unclear how the quality of current CL methods is
impacted by aggregation.

Figure 1: Overview of Freddie. An Orchestrator and an Extractor are
respectively used for deploying experiments and collecting data. Experiments
are run as TrainJobs managed by Kubeflow Training Operators. Within such a
job, the experiment is controlled by the federator and learned by the clients.

##  4 Freddie: A Framework for Reproducible FCL Research

In this section, for space reasons, we focus on Freddie’s implementation based
on containerization and orchestration methods, and on its support for FCL.

Kubernetes and containers. Fig. 1 presents a high-level overview of Freddie.
The Orchestrator provides functionality to kick-off experiments, in turn,
deployed on the cluster by Kubeflow’s training operators. The extractor
provides volumes for experiments to write to provided by an NFS provisioner
and server. Experiments themselves are performed by federator and client
nodes. The overall flow of the federated learning system with Freddie is as
follows. First, the user submits an experiment description of the system and
hyper-parameters. Following this, the Orchestrator deploys and monitors the
experiment. This design allows the user to scale his experiment up from small-
scale prototyping with minimal effort. The Extractor allows users to store and
retrieve experiment statistics and artefacts created by the federator or
clients.

The communication between any two parties in the system is asynchronous,
allowing the development of FL systems with non-blocking federator-client
interactions. With this flexibility, clients can run in Kubernetes clusters
and on individual machines. To allow users to describe their experiment as
distributed code, we rely on Kubeflow’s training operators, which provide a
means to set up distributed learning on Kubernetes with popular ML libraries.

Novel support for FCL. Freddie supports the SOTA algorithms for FCL [18] and
common CL methods such as EWC [5] and GEM [8]. For CL, Freddie implements
Task-IL and Domain-IL [16] through the use of sliding, expanding, and full
window mechanisms. A sliding window restricts the output classes only to those
of the task evaluated at a time $t$. Expanding windows do not utilize the task
IDs to make any such restriction, so the output classes include all classes
learned until that time. A full window does not restrict outputs based on the
task and can be used in the standard federated learning scenario.

The added complexity of FCL allows for workloads over the same set of tasks to
produce different results. We devise three different schemes that partition
tasks differently, and that can be used to evaluate a FCL scheme over a
representative set of scenarios. We discuss these three schemes, which we coin
Column, Balanced, and Shuffled, respectively. Column, as shown in Fig. 2(a),
all clients handle tasks in the same order. Thereby naively splitting the CL
workload across clients, resulting in a significant expected catastrophic
forgetting effect. Balanced, as depicted in Fig. 2(b), aims to prevent
catastrophic forgetting by organizing the tasks so that different clients
train on them across multiple steps. We propose a partitioning scheme to
lessen the effect of catastrophic forgetting, resulting in a task being
trained on by at most one client. While this scheme addresses short term-
forgetting, long-term forgetting may still occur. Shuffled, see Fig. 2(c),
randomly orders each client’s task, thereby relying on pseudo-randomness in
conjunction with a pre-specified seed.

(a) Column-wise.

(b) Balanced.

(c) Shuffled.

Figure 2: Column, balanced and shuffled task partition schemes for CL.

##  5 Performance evaluation

We demonstrate some features of Freddie through experiments. For this purpose,
we use the overlapping CIFAR100 dataset. The labels that already exist in
CIFAR100 are used to partition the data into different tasks for FCL,
following the same steps as in [18]. FL can be viewed as a corner case of FCL
where there is a single task consisting of the whole dataset. We first
consider a FL scenario with the default version of CIFAR100, and then consider
the overlapping CIFAR100 split into 10 separate tasks in a FCL scenario. We
use the average accuracy metric following the CL literature [2, 13].

Scalability. To investigate Freddie’s emulation capability, we perform a
_small_ and _large_ scale experiment on a Google Kubernetes Engine (GKE)
cluster to cover possible use cases. During deployment, the pods of the
federation and clients were run on a separate node pool scaled to meet each
experiment’s requirements. We study the performance of an FL experiment
emulated on a CPU-enabled Kubernetes cluster, where multiple clients may run
on a single Kubernetes node. Parameters of the experiments are provided in
Tab. 1.

Table 1: System and hyperparameters used in ‘small’ and ‘scale’ experiments. All experiments were run on ‘e2-standard-8’ nodes. | System | Federator (F) | Clients (C)
---|---|---|---
| Nodes | #C |  | CPU
---
(F/C)
| Memory
---
(F/C)
Strategy |  | #Rounds
---
(R)
#C/R | Model | Data | BS
Small | 2,2,3 | 5,10,20 | 2/1 | 2/2G | FedAvg | 100 | 5 | LeNet | CIFAR10 | 64
Scale | 4,12 | 25,75 | 2/2,6G | 85 | all | ResNet | CIFAR100

The _small_ experiment in Fig. 3(a) depicts the spread round times of clients
(scaled) and the federator, with 5 selected clients per round. The client
round duration is scaled by the number of clients (World Size WS)
($|\mathcal{D}_{Cifar}|/\text{WS}$) to account for differences in clients’
datasets as the WS increases. The outliers in the plots originate from the
first epoch run on clients, which are inherently slower due to loading data
into memory. Nevertheless, it is expected that the scaled client duration
stays relatively constant, while the result shows an increase as the number of
clients increases (from 115 to 123 s, and from 136 to 138 s). Similarly, the
federator sees a positive correlation between round duration and WS. The
number of co-scheduled clients on the same node can explain this trend, as the
networking overhead stays the same.

For the scale experiments, we provide the round time density estimate in Fig.
3(b). The client round times exhibit the same range of processing times that
were observed in the ‘small’ setting. In both settings, participating clients
in each round may run on the same node, varying from 3 to 7 clients per node.
We use similar settings in the ‘small’ configuration that involves 20 clients,
where 4 nodes are used. As such, confirming that resource contingencies due to
co-scheduling will likely cause the increased client round time with 20
clients. The different modes within the client’s round duration can be
explained by imperfect data splits and the imbalanced assignment of the number
of clients to be co-scheduled with the federator. The federator’s density
estimate shows a similar pattern with two distinct modes. With the cluster
configurations employed, i.e., 4 and 12 nodes, it is possible for the
federator to be co-scheduled on a machine with different numbers of clients.
As a result, the federator experiences a variable level of resource
contingency. However, an increase in the two modes is visible as the number of
clients increases, which is expected due to the increased communication
volumes.

(a) Experiments with 5 participating clients per round (n=1).

(b) Scaled experiments where all clients participate in each round (n=3).

Figure 3: Client and federator round durations with Freddie for small (5-20
clients, LeNet5 & CIFAR10) and large scale experiment (ResNet-18 & CIFAR100).
Client durations are scaled by the total number of clients (WS).

Task-IL vs Domain-IL. Assuming that the model knows the ID of the task it is
currently training on, or evaluating, increases its accuracy. For FCL, Task-
Interactive Learning and Domain-Interactive Learning are implemented using the
sliding and expanding-window, respectively. Let us recall that sliding-windows
use task IDs, contrary to expanding-windows. For the overlapping CIFAR100
dataset, if one assumes that the task ID is known, then the number of output
classes is restricted to only the five sub-classes within that task. Thus
leading to higher average task probabilities for Task-IL scenarios. This
difference is prevalent in Fig. 4(a). Under the expanding window scheme,
classification outputs one of $5T$ classes, where $T$ is the number of tasks
learned until evaluation time. Therefore, the probability of classifying
correctly is even lower than in the sliding window scenario. Fig. 4(a) shows
the positive impact of using a task ID on accuracy. Using sliding-window
results in higher accuracy than expanding-window, which sometimes has to be
used because of the application use case. Because of this difference, Freddie
supports both Task-IL and Domain-IL.

(a) Reduced catastrophic forgetting effect window (Task-IL) compared to
expanding window (Domain-IL).

(b) Impact of partition scheme on evaluation accuracy (10 tasks of overlapping
CIFAR100, LeNet-5, 10 rounds of 2 epochs per task).

Figure 4: Impact of task heterogeneity on FCL.

FCL Task Heterogeneity. As discussed in Section 3, tasks can be processed in
different orders at each client. To demonstrate the different effects that
different sequences of tasks produce, we implement the Overlapped-CIFAR100
dataset with 20 tasks that can be used for FCL [18]. The accuracy in Fig. 4(b)
is calculated as the average accuracy of all tasks seen until that point,
resulting in expected ‘drops’ in accuracy as new tasks are introduced. Indeed,
the learning curves in Fig. 4(b) show noticeable drops over time. However,
different trends are visible between workloads. The column scheme suffers more
from more pronounced _catastrophic forgetting_ than the shuffled and balanced
scheme, resulting in lower accuracy. We observe that the column scheme, on
average, results in a 4% test accuracy drop compared to the column and
shuffled schemes.

##  6 Conclusion

We presented Freddie, the first framework for reproducible Federated Continual
Learning research, which is motivated by the increasing importance of
Federated and Continual Learning. Freddie’s deployment abilities on different
platforms, scalability with the number of clients, and support for data and
task heterogeneity provide FL practitioners with a powerful tool. Our
experimental results showcase previously unaddressed performance issues that
Federated Continual Learning systems might face: severe catastrophic
forgetting in different task heterogeneity settings. Freddie is open-source,
and will soon be extended to support new CL datasets, algorithms, and
generative models.

## References

  * [1] Beutel, D.J., Topal, T., Mathur, A., et al.: Flower: A friendly federated learning research framework. CoRR abs/2007.14390 (2020)
  * [2] Chaudhry, A., Ranzato, M., Rohrbach, M., et al.: Efficient lifelong learning with A-GEM. In: ICLR (2019)
  * [3] FedAI: Fate (2019), https://fate.fedai.org/
  * [4] Ke, Z., Liu, B., Ma, N., et al.: Achieving forgetting prevention and knowledge transfer in continual learning. In: NeurIPS (2021)
  * [5] Kirkpatrick, J., Pascanu, R., Rabinowitz, N.C., et al.: Overcoming catastrophic forgetting in neural networks. CoRR abs/1612.00796 (2016)
  * [6] Lange, M.D., Aljundi, R., Masana, M., et al.: Continual learning: A comparative study on how to defy forgetting in classification tasks. CoRR abs/1909.08383 (2019)
  * [7] Lomonaco, V., Pellegrini, L., Cossu, A., et al.: Avalanche: an end-to-end library for continual learning. In: CVPR (2021)
  * [8] Lopez-Paz, D., Ranzato, M.A.: Gradient episodic memory for continual learning. In: NeurIPS (2017)
  * [9] Ma, Y., Yu, D., Wu, T., et al.: PaddlePaddle: An open-source deep learning platform from industrial practice. Frontiers of Data and Computing 1(1), 105–115 (2019)
  * [10] Mallya, A., Lazebnik, S.: Packnet: Adding multiple tasks to a single network by iterative pruning. In: CVPR (2018)
  * [11] Masana, M., Liu, X., Twardowski, B., Menta, M., et al.: Class-incremental learning: Survey and performance evaluation on image classification. IEEE Trans. Pattern Anal. Mach. Intell. 45(5), 5513–5533 (2023)
  * [12] McMahan, B., Moore, E., Ramage, D., et al.: Communication-efficient learning of deep networks from decentralized data. In: Artificial intelligence and statistics. pp. 1273–1282. PMLR (2017)
  * [13] Mirzadeh, S.I., Farajtabar, M., Pascanu, R., et al.: Understanding the role of training regimes in continual learning. NeurIPS (2020)
  * [14] Reina, G.A., Gruzdev, A., Foley, P., et al.: Openfl: An open-source framework for federated learning. CoRR abs/2105.06413 (2021)
  * [15] Shin, H., Lee, J.K., Kim, J., Kim, J.: Continual learning with deep generative replay. In: NeurIPS. pp. 2990–2999 (2017)
  * [16] van de Ven, G.M., Tolias, A.S.: Three scenarios for continual learning. CoRR abs/1904.07734 (2019)
  * [17] Wołczyk, M., Zajkac, M., Pascanu, R., et al.: Continual world: A robotic benchmark for continual reinforcement learning. In: NeurIPS (2021)
  * [18] Yoon, J., Jeong, W., Lee, G., et al.: Federated continual learning with weighted inter-client transfer. In: ICML (2021)
  * [19] Zhou, D.W., Wang, F.Y., Ye, H.J., et al.: Pycil: A python toolbox for class-incremental learning. CoRR abs/2112.12533 (2021)


# The Fictitious Domain Method Based on Navier Slip Boundary Condition for
Simulation of Flow-Particle Interaction

Rong Zhang  Qiaolin He  Key Lab of Intelligent Analysis and Decision on
Complex Systems, Chongqing University of Posts and Telecommunications,
Chongqing 400065, China<EMAIL_ADDRESS> School of Mathematics,
Sichuan University, Chengdu, 610064, China<EMAIL_ADDRESS>
###### Abstract

In this article, we develop a least–squares/fictitious domain method for
direct simulation of fluid particle motion with Navier slip boundary condition
at the fluid–particle interface. Let $\Omega$ and $B$ be two bounded domains
of $\mathbb{R}^{d}$ such that $\overline{B}\subset\Omega$. The motion of solid
particle $B$ is governed by Newton’s equations. Our goal here is to develop a
fictitious domain method where one solves a variant of the original problem on
the full $\Omega$, followed by a well–chosen correction over $B$ and
corrections related to translation velocity and angular velocity of the
particle. This method is of the virtual control type and relies on a
least–squares formulation making the problem solvable by a conjugate gradient
algorithm operating in a well chosen control space. Since the fully explicit
scheme to update the particle motion using Newton’s equation is unstable, we
propose and implement an explicit–implicit scheme in which, at each time step,
the position of the particle is updated explicitly, and the solution of
Navier-Stokes equations and particle velocities are solved by the the
least–squares/fictitious domain method implicitly. Numerical results are given
to verify our numerical method.

###### keywords:

least–squares, fictitious domain method, incompressible viscous flow, Navier
slip boundary condition

††journal: Elsevier

##  1 Introduction

In order to understand the interactions between the solid particle and the
fluid, we take a more fundamental approach with the Navier–Stokes equations
for fluid and Newton’s equations of motion for the particle. Numerical
simulation of these equations can give us all the details of the flow and the
fluid–particle interactions together with clear understandings of the
mechanisms involved. There are some numerical works on these types of
simulations. Authors in [13] and [11] used arbitrary Lagrangian–Eulerian
method to directly simulate the motion of fluid and particles, which based on
moving unstructured grids and needs remeshing and projection.

Another important method is fictitious domain method. In [5, 7, 6, 8],
fictitious domain methods were discussed for the solution of Dirichlet
problems, where the Dirichlet boundary condition being enforced as a side
constraint, using a boundary supported Lagrange multiplier. A volume–supported
Lagrange multiplier based fictitious domain method was introduced in [14], the
main motivation being the direct numerical simulation of particulate flow when
the number of particles exceeds $10^{3}$. Initially tested on particulate flow
with spherical particles, the method discussed in [14] was generalized to
situations involving particles with more complicated shapes, as shown for
example in [15]. A (brief) history of fictitious domain methods  can be found
in, e.g., Chapter 8 in [2]. Examples of non Lagrange multiplier based
fictitious domain methods can be found in the immersed boundary method for the
simulation of incompressible viscous flow in regions with elastic moving
boundaries [16]. Then, M. Uhlmann [17] used the immersed boundary method to
simulate particulate flow.

The main idea behind fictitious domain methods is to extend a problem
initially posed on a geometrically complex shaped domain to a larger simpler
domain; this provides two main advantages when constructing numerical schemes:
(i) The extended domain is geometrically simpler and allows the use of fast
solvers. (ii) The same fixed mesh can be used for the entire computation,
eliminating thus the need for repeated remeshing and projection. All the
studies that we know, concerning the application of fictitious domain methods
to the simulation of particulate flow, consider no-slip boundary conditions at
the interface between fluid and particles. There are situations however, in
micro-fluidics for example, where a slip condition on the particle surface is
more realistic than the no-slip one. If the no-slip boundary condition on the
particle surface is replaced by the Navier slip boundary condition, the
volume–supported Lagrange multiplier based fictitious domain methods discussed
in [5, 6, 7, 8, 14, 15], which rely on $H^{1}$– extensions, are not easy to
generalize to the slip situation. For incompressible viscous flow around
obstacle with Navier slip boundary condition [10], we have developed a new
least–squares/fictitious domain method. The method in [10] is of the virtual
control type (in the sense of J. L. Lions; see [3]) and relies on a
least–squares formulation making the problem solvable by a conjugate gradient
algorithm operating in a well–chosen control space, which is a generation of
least–squares/fictitious domain method for linear elliptic problem with Robin
boundary condition (see [4]) to Navier–Stokes problem with Navier slip
boundary condition. Such least–squares idea has been sucessfully generated to
Navier–Stokes–Cahn–Hilliard system in [9].

The main goal of the present article is to discuss the numerical simulation of
fluid-particle interaction when Navier slip take effect at the interface
fluid/particle. The particle velocity is updated together with the solution of
Navier–Stokes equations for the stability. The position of particle is updated
explicitly. The formulation of the problem is given in Section 2. In Section
3, we describe an operator–splitting scheme for the problem, which including
the discription of a least–squares/fictitious domain method for the sub-
problem and a conjugate gradint solution of the least–squares problem. The
finite element implementation of the above methodology is discussed in Section
4. Finally, we present in Section 5 the results of numerical experiments.

##  2 Problem formulation

Let $\Omega$ be a bounded domain of $\mathbb{R}^{d}(d=2$ or 3), including the
interior of the particle. We denote by $\Gamma$ the boundary of $\Omega$. The
geometry of the problem has been visualized on Figure 1. We suppose that
$\Omega$ contains a rigid body $B$ whose boundary is denoted by $\gamma$. For
simplicity, we assume that the fluid velocity satisfies a Dirichlet boundary
condition on the outer boundary $\Gamma$. The $\rho_{s}$ and $\rho_{f}$ are
particle density and fluid density, respectively. Compared to previous work,
the main novelty here is that we will assume that a slip boundary condition
holds at the interface $\gamma$ between $\Omega$ and $B$. Assuming that the
external force is $\mathbf{g}$, the motion of fluid–particle mixture is
governed by the following equations:
Fluid motion

$\displaystyle\rho_{f}\left[\frac{\partial\mathbf{u}}{\partial
t}+(\mathbf{u}\cdot\nabla)\mathbf{u}\right]=\rho_{f}\mathbf{g}+\nabla\cdot\boldsymbol{\sigma},\quad\mbox{in}\quad\Omega\setminus\overline{B(t)}$
(2.1)
$\displaystyle\nabla\cdot\mathbf{u}=0\quad\mbox{in}\quad\Omega\setminus\overline{B(t)}$
(2.2)
$\displaystyle(\mathbf{u}-\mathbf{u}_{B})\cdot\mathbf{n}=0\quad\mbox{on}\quad\gamma\times(0,T),$
(2.3)
$\displaystyle(\boldsymbol{\sigma}\mathbf{n}+\frac{\mu}{l_{s}}(\mathbf{u}-\mathbf{u}_{B}))\times\mathbf{n}=0\quad\mbox{on}\quad\gamma\times(0,T),$
(2.4)
$\displaystyle\mathbf{u}=\mathbf{g}_{\Gamma}\quad\mbox{on}\quad\Gamma\times(0,T),$
(2.5)
$\displaystyle\mathbf{u}_{B}=\mathbf{U}+\boldsymbol{\omega}\times\mathbf{r}\
\mbox{on}\ \gamma,$ (2.6)
$\displaystyle\mathbf{u}(x,0)=\mathbf{u}_{0},x\in\Omega\setminus\overline{B(0)},\quad\mbox{with}\quad\nabla\cdot\mathbf{u_{0}}=0.$
(2.7)

Particle motion

$\displaystyle\mathbf{M}\frac{d\mathbf{U}}{dt}=\mathbf{M}\mathbf{g}+\mathbf{F},$
(2.8) $\displaystyle\frac{d\mathbf{I}\boldsymbol{\omega}}{dt}=\mathbf{T},$
(2.9) $\displaystyle\mathbf{U}|_{t=0}=\mathbf{U}_{0},$ (2.10)
$\displaystyle\boldsymbol{\omega}|_{t=0}=\boldsymbol{\omega}_{0}$ (2.11)
$\displaystyle\frac{d\mathbf{X}}{dt}=\mathbf{U},$ (2.12)
$\displaystyle\mathbf{X}|_{t=0}=\mathbf{X}_{0},$ (2.13)

where $\mathbf{n}$ denotes the outward normal unit vector at
$\partial(\Omega\backslash\overline{B})$, pointing outward to the fluid
region, $\mathbf{r}\triangleq\mathbf{x}-\mathbf{X}$; $\mathbf{u}$ and
$\boldsymbol{\sigma}$ are the fluid velocity and stress;
$\mathbf{M},\mathbf{I},\mathbf{U},\boldsymbol{\omega},\mathbf{X}$ are the
mass, moment of inertia,translational velocity, angular velocity, center of
mass of particle respectively; and $\mathbf{F}=-\int_{\partial
B(t)}\boldsymbol{\sigma}\mathbf{n}d\gamma$, $\mathbf{T}=-\int_{\partial
B(t)}\mathbf{r}\times\boldsymbol{\sigma}\mathbf{n}d\gamma$ are the
hydrodynamics force and torque on the particle, and $l_{s}>0$ is the slip
length, $\mathbf{u}_{B}$ is the velocity of the rigid body, $\mu$ is positive
constant and denote fluid viscosity. Here
$\boldsymbol{\sigma}=2\mu\mathbf{D}(\mathbf{u})-p\mathbf{E}$ with
$\mathbf{D}(\mathbf{v})=\frac{1}{2}[\nabla\mathbf{v}+(\nabla\mathbf{v})^{t}],\forall\
\mathbf{v}$, where $\mathbf{E}$ is identity matrix.

A classical variational formulation of the problem (2.1)–(2.7) is given by:
Find $\\{\mathbf{u},p,\mathbf{U},\boldsymbol{\omega}\\}$ such that, a.e. on
$(0,T)$, (2.14)–(2.18) are established.

$\displaystyle\rho_{f}\int_{\Omega\backslash\bar{B}}\left[\frac{\partial\mathbf{u}}{\partial
t}+(\mathbf{u}\cdot\nabla)\mathbf{u}\right]\cdot\mathbf{v}dx+2\mu\int_{\Omega\backslash\bar{B}}\mathbf{D}(\mathbf{u}):\mathbf{D}(\mathbf{v})dx$
$\displaystyle+\int_{\Omega\backslash\bar{B}}\nabla
p\cdot\mathbf{v}dx+\frac{\mu}{l_{s}}\int_{\gamma}(\mathbf{u}-(\mathbf{U}+\boldsymbol{\omega}\times\mathbf{r}))\cdot\mathbf{v}d\gamma=\rho_{f}\int_{\Omega\backslash\bar{B}}\mathbf{g}\cdot\mathbf{v}dx,$
(2.14) $\displaystyle\forall\mathbf{v}\in\mathbb{V}_{0},$
$\displaystyle\int_{\Omega\backslash\bar{B}}\nabla\cdot\mathbf{u}\ qdx=0,\
\forall q\in L^{2}(\Omega\backslash\bar{B}),$ (2.15)
$\displaystyle\mathbf{u}=\mathbf{g}_{\Gamma}\ \mbox{on}\ \Gamma,$ (2.16)
$\displaystyle\mathbf{u}(0)=\mathbf{u}_{0},$ (2.17)
$\displaystyle\left(\mathbf{u}-(\mathbf{U}+\boldsymbol{\omega}\times\mathbf{r})\right)\cdot\mathbf{n}=0\
\mbox{on}\ \gamma\times(0,T),$ (2.18)

where
$\mathbb{V}_{0}=\\{\mathbf{v}|\mathbf{v}\in(H^{1}(\Omega\backslash\bar{B}))^{d},\
\mathbf{v}=\mathbf{0}\ \mbox{on}\ \Gamma,\ \mathbf{v}\cdot\mathbf{n}=0\
\mbox{on}\ \gamma\\}.$

Figure 1: Problem geometry

##  3 An operator–splitting scheme for the time–discretization of problem
(2.14)–(2.18) and (2.8)–(2.13)

We first came up with a very simple scheme which decoupled the motion of the
particle and the motion of fluid at each time step. That is, we use the method
in [10] to solve (3.1)–(3.4) firstly, solve the pure advection problem
(3.5)–(3.8) secondly, update particle velocity and position explicitly.
Unfortunately, this simple scheme is unstable. For $n\geq 0$, $\mathbf{u}^{n}$
being known, $\Delta t>0$ is a time–discretization step and $t^{n}=n\Delta t$,

$\displaystyle\rho_{f}\int_{\Omega\backslash\bar{B}}\left[\frac{\mathbf{u}^{n+\frac{1}{2}}-\mathbf{u}^{n}}{\Delta
t}\cdot\mathbf{v}dx\right]+2\mu\int_{\Omega\backslash\bar{B}}\mathbf{D}(\mathbf{u}^{n+\frac{1}{2}}):\mathbf{D}(\mathbf{v})dx+\int_{\Omega\backslash\bar{B}}\nabla
p^{n+1}\cdot\mathbf{v}dx$
$\displaystyle+\frac{\mu}{l_{s}}\int_{\gamma}(\mathbf{u}^{n+\frac{1}{2}}-\mathbf{u}_{B})\cdot\mathbf{v}d\gamma=\rho_{f}\int_{\Omega\backslash\bar{B}}\mathbf{g}\cdot\mathbf{v}dx,\
\forall\mathbf{v}\in\mathbb{V}_{0},$ (3.1)
$\displaystyle\int_{\Omega\backslash\bar{B}}\nabla\cdot\mathbf{u}^{n+\frac{1}{2}}qdx=0,\
\forall q\in L^{2}(\Omega\backslash\bar{B}),$ (3.2)
$\displaystyle(\mathbf{u}^{n+\frac{1}{2}}-\mathbf{u}_{B})\cdot\mathbf{n}=0\
\mbox{on}\ \gamma,$ (3.3)
$\displaystyle\mathbf{u}^{n+\frac{1}{2}}=\mathbf{g}_{\Gamma}\ \mbox{on}\
\Gamma,$ (3.4)

where $\mathbf{u}_{B}=\mathbf{U}+\boldsymbol{\omega}\times\mathbf{r}$ is
obtained by (2.8) and (2.9). In order to solve problem (2.14)–(2.18) and
(2.8)–(2.13) numerically, we advocate the operator–splitting scheme where
fluid velocity and particle velocty are updated implicitly together and
particle position is updated explicitly. Then we solve the following pure
advection problem

$\displaystyle\frac{\partial\mathbf{u}}{\partial
t}+(\mathbf{u}^{n+\frac{1}{2}}\cdot\nabla)\mathbf{u}=\mathbf{0}\ \mbox{in}\
(\Omega\backslash\bar{B})\times(t^{n},t^{n+1}),$ (3.5)
$\displaystyle\mathbf{u}(t^{n})=\mathbf{u}^{n+\frac{1}{2}},$ (3.6)
$\displaystyle\mathbf{u}=\mathbf{g}_{\Gamma}\ \mbox{on}\
\mathbf{\Gamma}_{-}\times(t^{n},t^{n+1});$ (3.7)
$\displaystyle\mathbf{u}=\mathbf{u}^{n+\frac{1}{2}}\ \mbox{on}\
\mathbf{\gamma}_{-}\times(t^{n},t^{n+1}),$ (3.8)

set

$\displaystyle\mathbf{u}^{n+1}=\mathbf{u}(t^{n+1}).$ (3.9)

where,
$\mathbf{\Gamma}_{-}=\\{x|x\in\Gamma,\mathbf{g}_{\Gamma}\cdot\mathbf{n}<0\\}$
and $\mathbf{\gamma}_{-}=\\{x|x\in\gamma,\mathbf{u}_{B}\cdot\mathbf{n}<0\\}$.
Since we consider low Reynolds number flows, the backward method of
characteristics will be used to solve problems (3.5)–(3.8). The implementation
of such method has been detailed described in Chapter 6 in [2]. Finally, the
position of particle is updated by the approach in [15] for (2.12)–(2.13). At
each time step, keeping the distance constant between points $\mathbf{X}_{i}$
in particle is important since we are dealing with rigid particle. Therefore,
our paper will focus on how to solve (3.1)–(3.4) and updating
$\mathbf{u}_{B}$.

###  3.1 On the fictitious domain solution of sub-problem excluding advection
terms

In order to solve fluid velocity and particle velocity together, we advocate
the following variant of _virtual control/ fictitious domain method_ discussed
in [4] and [10] (below, $\mathbf{u}_{\ast},\ p_{1\ast}$ are defined over the
whole $\Omega$ and verify
$\mathbf{u}_{\ast}|_{\Omega\backslash\bar{B}}=\mathbf{u}^{n},\
p_{1\ast}|_{\Omega\backslash\bar{B}}=p^{n},\mathbf{U}_{\ast}=\mathbf{U}^{n},\omega_{\ast}=\omega^{n}$
). For simlicity, we consider two dimensional motion (the following numerical
method can be generated to three dimension case straightforwardly), therefore
angular velocity $\boldsymbol{\omega}$ reduced to $\omega$ and angular
velocity equation becomes

$\displaystyle\mathbf{I}\frac{d\omega}{dt}=\mathbf{T}.$

Suppose there exists
$\mathbf{y}\in(L^{2}(B))^{2},\mathbf{C}_{1}\in\mathbb{R}^{2},C_{2}\in\mathbb{R}$
such that the following relations hold:

$\displaystyle\mathbf{u}_{1}\in(H^{1}(\Omega))^{2},\mathbf{u}_{1}=\mathbf{g}_{\Gamma}\
\mbox{on}\ \Gamma,\ p_{1}\in H^{1}(\Omega),$
$\displaystyle\rho_{f}\int_{\Omega}\frac{\mathbf{u}_{1}-\mathbf{u}_{\ast}}{\Delta
t}\cdot\mathbf{v}dx-\rho_{f}\int_{B}\frac{\mathbf{y}}{\Delta
t}\cdot\mathbf{v}dx+2\mu\int_{\Omega}\mathbf{D}(\mathbf{u}_{1}):\mathbf{D}(\mathbf{v})dx+\int_{\Omega}\nabla
p_{1}\cdot\mathbf{v}dx$
$\displaystyle=\rho_{f}\int_{\Omega\backslash\bar{B}}\mathbf{g}\cdot\mathbf{v}dx,\
\forall\mathbf{v}\in(H_{0}^{1}(\Omega))^{2},$ (3.10)
$\displaystyle\int_{\Omega}\nabla\cdot\mathbf{u}_{1}qdx=-\int_{\Omega}\mathbf{u}_{1}\cdot\nabla
qdx+\int_{\Gamma}(\mathbf{g}_{\Gamma}\cdot\mathbf{n})qd\Gamma=0,\ \forall q\in
H^{1}(\Omega),$ (3.11) $\displaystyle\mathbf{u}_{2}\in(H^{1}(B))^{2},\
p_{2}\in H^{1}(B),\
\left(\mathbf{u}_{2}-(\widetilde{\mathbf{U}}+\widetilde{\omega}\times\mathbf{r})-\frac{1}{\Delta
t}\left(\mathbf{C}_{1}+C_{2}\times\mathbf{r}\right)\right)\cdot\mathbf{n}=0\
\mbox{on}\ \gamma,$
$\displaystyle\rho_{f}\int_{B}\frac{\mathbf{u}_{2}-\mathbf{u}_{\ast}-\mathbf{y}}{\Delta
t}\cdot\mathbf{v}dx+2\mu\int_{B}\mathbf{D}(\mathbf{u}_{2}):\mathbf{D}(\mathbf{v})dx+\int_{B}\nabla
p_{2}\cdot\mathbf{v}dx$
$\displaystyle=\frac{\mu}{l_{s}}\int_{\gamma}(\mathbf{u}_{1}-(\widetilde{\mathbf{U}}+\widetilde{\omega}\times\mathbf{r}))\cdot\mathbf{v}d\gamma-\frac{1}{\Delta
t}\frac{\mu}{l_{s}}\int_{\gamma}(\mathbf{C}_{1}+C_{2}\times\mathbf{r})\cdot\mathbf{v}d\gamma,$
(3.12)
$\displaystyle\forall\mathbf{v}\in\mathbb{V}_{0B}\triangleq\\{(H^{1}(B))^{2},\mathbf{v}\cdot\mathbf{n}=0\
\mbox{on}\ \gamma\\},$
$\displaystyle\int_{B}\nabla\cdot\mathbf{u}_{2}qdx=-\int_{B}\mathbf{u}_{2}\cdot\nabla
qdx+\int_{\gamma}(\widetilde{\mathbf{U}}+\widetilde{\omega}\times\mathbf{r})\cdot\mathbf{n}qd\gamma+\frac{1}{\Delta
t}\int_{\gamma}(\mathbf{C}_{1}+C_{2}\times\mathbf{r})\cdot\mathbf{n}qd\gamma$
$\displaystyle=0,\ \forall q\in H^{1}(B),$ (3.13)
$\displaystyle\mathbf{U}=\widetilde{\mathbf{U}}+\frac{1}{\Delta
t}\mathbf{C}_{1},$ (3.14)
$\displaystyle\omega=\widetilde{\omega}+\frac{1}{\Delta t}C_{2},$ (3.15)

where

$\displaystyle\widetilde{\mathbf{U}}=\mathbf{U}_{\ast}-\frac{\Delta
t}{\mathbf{M}}\int_{\gamma}\left(-p_{1\ast}\mathbf{E}+2\mu\mathbf{D}(\mathbf{u}_{\ast})\right)\mathbf{n}d\gamma,$

and

$\displaystyle\widetilde{\omega}=\omega_{\ast}-\frac{\Delta
t}{\mathbf{I}}\int_{\gamma}\mathbf{r}\times\left(-p_{1\ast}\mathbf{E}+2\mu\mathbf{D}(\mathbf{u}_{\ast})\right)\mathbf{n}d\gamma.$

Both problems (3.10)–(3.11) and (3.12)–(3.13) have a unique solution in
$(H^{1}(\Omega))^{2}\times H^{1}(\Omega)/\mathbb{R}$ and $(H^{1}(B))^{2}\times
H^{1}(B)/\mathbb{R}$, respectively. We define
$\mathbf{A}:(L^{2}(B))^{2}\times\mathbb{R}^{2}\times\mathbb{R}\rightarrow(H^{1}(B))^{2}$
by

$\displaystyle\mathbf{A}(\mathbf{y},\mathbf{C}_{1},C_{2})=\left(\mathbf{u}_{2}-\mathbf{u}_{1}\right)\mid_{B}.$
(3.16)

Operator $\mathbf{A}$ is clearly affine and continuous. We observe that if
$\mathbf{y}^{\ast},\mathbf{C}^{\ast}_{1},C^{\ast}_{2}$ verify
$\mathbf{A}(\mathbf{y},\mathbf{C}_{1},C_{2})=\mathbf{0}$, we then have
$\mathbf{u}_{2}=\mathbf{u}_{1}$ on $B$ and it is easy to see that
$\mathbf{u}_{1}|_{\Omega\setminus\bar{B}},p_{1}|_{\Omega\setminus\bar{B}},\mathbf{U}$
and $\omega$ are what we need. We will numerically discuss the solution of the
following functional equation

$\displaystyle\mathbf{A}(\mathbf{y}^{\ast},\mathbf{C}^{\ast}_{1},C^{\ast}_{2})=\mathbf{0}.$
(3.17)

######  Remark 3.1.

Problem (3.17) can be viewed as an exact controllability problem in the sense
of [3].

####  3.1.1  A least–squares formulation

In order to solve system (3.17), we use the following least–squares approach:

$\displaystyle\mbox{Find}\
\mathbf{y}^{\ast}\in(L^{2}(B))^{2},\mathbf{C}_{1}^{\ast},C_{2}^{\ast}\
\mbox{such that}$ $\displaystyle
J(\mathbf{y}^{\ast},\mathbf{C}^{\ast}_{1},C^{\ast}_{2})\leq
J(\mathbf{y},\mathbf{C}_{1},C_{2}),\
\forall\mathbf{y}\in(L^{2}(B))^{2},\mathbf{C}_{1}\in\mathbb{R}^{2},\
C_{2}\in\mathbb{R},$ (3.18)

where

$\displaystyle
J(\mathbf{y},\mathbf{C}_{1},C_{2})=\frac{1}{2}\left[\rho_{f}\int_{B}|\mathbf{u}_{2}-\mathbf{u}_{1}|^{2}dx+2\mu\Delta
t\int_{B}|\mathbf{D}(\mathbf{u}_{2}-\mathbf{u}_{1})|^{2}dx\right],$ (3.19)

where in (3.19), $\mathbf{u}_{1}$ and $\mathbf{u}_{2}$ are obtained from
$\mathbf{y},\mathbf{C}_{1}$ and $C_{2}$, via the solutions of (3.10)–(3.11)
and (3.12)–(3.13). A natural candidate for the solution of the minimization
problem (3.18)–(3.19) is a conjugate gradient algorithm operating in the
virtual control space $(L^{2}(B))^{2}\times\mathbb{R}^{2}\times\mathbb{R}$.
Such a solution is characterized by

$\displaystyle\frac{\partial
J}{\partial\mathbf{y}}(\mathbf{y}^{\ast},\mathbf{C}^{\ast}_{1},C^{\ast}_{2})=\mathbf{0},\
\frac{\partial
J}{\partial\mathbf{C}_{1}}(\mathbf{y}^{\ast},\mathbf{C}^{\ast}_{1},C^{\ast}_{2})=\mathbf{0},\
\frac{\partial J}{\partial
C_{2}}(\mathbf{y}^{\ast},\mathbf{C}^{\ast}_{1},C^{\ast}_{2})=0.$

####  3.1.2 On the computation of partial differentials of functional $J$

We are going to address this most important issue using a _perturbation
approach_ (as done in [4] for a close variant of problem (3.18), and in [3]
for various linear and nonlinear control problems). Suppose that
$\mathbf{y}\in(L^{2}(B))^{2},\mathbf{C}_{1}\in\mathbb{R}^{2}$ and
$C_{2}\in\mathbb{R}$, a _perturbation_ $\delta\mathbf{y}$ also belonging to
$(L^{2}(B))^{2},\delta\mathbf{C}_{1}$ belonging to $\mathbb{R}^{2}$ and
$\delta C_{2}$ belonging to $\mathbb{R}$ are given, we have then

$\displaystyle\delta
J(\mathbf{y},\mathbf{C}_{1},C_{2})=\rho_{f}\int_{B}(\mathbf{u}_{2}-\mathbf{u}_{1})\cdot\delta(\mathbf{u}_{2}-\mathbf{u}_{1})dx+2\mu\Delta
t\int_{B}\mathbf{D}(\mathbf{u}_{2}-\mathbf{u}_{1}):\mathbf{D}\delta(\mathbf{u}_{2}-\mathbf{u}_{1})dx$
$\displaystyle=\rho_{f}\int_{B}(\mathbf{u}_{2}-\mathbf{u}_{1})\delta\mathbf{u}_{2}dx+2\mu\Delta
t\int_{B}\mathbf{D}(\mathbf{u}_{2}-\mathbf{u}_{1}):\mathbf{D}(\delta\mathbf{u}_{2})dx$
$\displaystyle+\rho_{f}\int_{B}(\mathbf{u}_{1}-\mathbf{u}_{2})\delta\mathbf{u}_{1}dx+2\mu\Delta
t\int_{B}\mathbf{D}(\mathbf{u}_{1}-\mathbf{u}_{2}):\mathbf{D}(\delta\mathbf{u}_{1})dx,$
(3.20)

with $\delta\mathbf{u}_{1}$ and $\delta\mathbf{u}_{2}$ verifying

$\displaystyle\delta\mathbf{u}_{1}\in(H^{1}_{0}(\Omega))^{2},\ \delta p_{1}\in
H^{1}(\Omega),$
$\displaystyle\rho_{f}\int_{\Omega}\delta\mathbf{u}_{1}\cdot\mathbf{v}dx+2\mu\Delta
t\int_{\Omega}\mathbf{D}(\delta\mathbf{u}_{1}):\mathbf{D}(\mathbf{v})dx+\Delta
t\int_{\Omega}\nabla\delta p_{1}\cdot\mathbf{v}dx$
$\displaystyle=\rho_{f}\int_{B}\delta\mathbf{y}\cdot\mathbf{v}dx,\
\forall\mathbf{v}\in(H^{1}_{0}(\Omega))^{2},$ (3.21)
$\displaystyle\int_{\Omega}\delta\mathbf{u}_{1}\cdot\nabla qdx=0,\ \forall
q\in H^{1}(\Omega),$ (3.22)
$\displaystyle\delta\mathbf{u}_{2}\in(H^{1}(B))^{2},\ \delta p_{2}\in
H^{1}(B),\ \left(\delta\mathbf{u}_{2}-\frac{1}{\Delta
t}\left(\delta\mathbf{C}_{1}+\delta
C_{2}\times\mathbf{r}\right)\right)\cdot\mathbf{n}=0\ \mbox{on}\ \gamma,$
$\displaystyle\rho_{f}\int_{B}\delta\mathbf{u}_{2}\cdot\mathbf{v}dx+2\mu\Delta
t\int_{B}\mathbf{D}(\delta\mathbf{u}_{2}):\mathbf{D}(\mathbf{v})dx+\Delta
t\int_{B}\nabla\delta p_{2}\cdot\mathbf{v}dx$ $\displaystyle=\frac{\mu\Delta
t}{l_{s}}\int_{\gamma}\delta\mathbf{u}_{1}\cdot\mathbf{v}d\gamma-\frac{\mu}{l_{s}}\int_{\gamma}\left(\delta\mathbf{C}_{1}+\delta
C_{2}\times\mathbf{r}\right)\cdot\mathbf{v}d\gamma+\rho_{f}\int_{B}\delta\mathbf{y}\cdot\mathbf{v}dx,\
\forall\mathbf{v}\in\mathbb{V}_{0B},$ (3.23)
$\displaystyle\int_{B}\delta\mathbf{u}_{2}\cdot\nabla qdx=\frac{1}{\Delta
t}\int_{\gamma}(\delta\mathbf{C}_{1}+\delta
C_{2}\times\mathbf{r})\cdot\mathbf{n}qd\gamma,\ \forall q\in H^{1}(B).$ (3.24)

Let us define the function $\boldsymbol{\pi}_{2},\ p_{\boldsymbol{\pi}_{2}}$
as the solution of the following linear variational problem:

$\displaystyle\boldsymbol{\pi}_{2}\in\mathbb{V}_{0B},\
p_{\boldsymbol{\pi}_{2}}\in H^{1}(B),$
$\displaystyle\rho_{f}\int_{B}\boldsymbol{\pi}_{2}\cdot\mathbf{v}dx+2\mu\Delta
t\int_{B}\mathbf{D}(\boldsymbol{\pi}_{2}):\mathbf{D}(\mathbf{v})dx+\Delta
t\int_{B}\nabla p_{\boldsymbol{\pi}_{2}}\cdot\mathbf{v}dx$
$\displaystyle=\rho_{f}\int_{B}(\mathbf{u}_{2}-\mathbf{u}_{1})\cdot\mathbf{v}dx+2\mu\Delta
t\int_{B}\mathbf{D}(\mathbf{u}_{2}-\mathbf{u}_{1})\cdot\mathbf{D}(\mathbf{v})dx,\
\forall\mathbf{v}\in\mathbb{V}_{0B},$ (3.25)
$\displaystyle\int_{B}\boldsymbol{\pi}_{2}\cdot\nabla qdx=0,\ \forall q\in
H^{1}(B).$ (3.26)

Taking $\mathbf{v}=\delta\mathbf{u}_{2}-\frac{1}{\Delta
t}\left(\delta\mathbf{C}_{1}+\delta C_{2}\times\mathbf{r}\right)$ in (3.25)
and $\mathbf{v}=\boldsymbol{\pi}_{2}$ in (3.23), using
$\mathbf{D}(\delta\mathbf{C}_{1}+\delta C_{2}\times\mathbf{r})=\mathbf{0}$,
and combining with (3.20), we obtain

$\displaystyle\delta
J(\mathbf{y},\mathbf{C}_{1},C_{2})=\rho_{f}\int_{B}\boldsymbol{\pi}_{2}\cdot\left(\delta\mathbf{u}_{2}-\frac{1}{\Delta
t}(\delta\mathbf{C}_{1}+\delta C_{2}\times\mathbf{r})\right)dx$
$\displaystyle+2\mu\Delta
t\int_{B}\mathbf{D}(\boldsymbol{\pi}_{2}):\mathbf{D}\left(\delta\mathbf{u}_{2}-\frac{1}{\Delta
t}(\delta\mathbf{C}_{1}+\delta C_{2}\times\mathbf{r})\right)dx$
$\displaystyle+\Delta t\int_{B}\nabla
p_{\boldsymbol{\pi}_{2}}\cdot\left(\delta\mathbf{u}_{2}-\frac{1}{\Delta
t}(\delta\mathbf{C}_{1}+\delta
C_{2}\times\mathbf{r})\right)dx+\frac{\rho_{f}}{\Delta
t}\int_{B}(\mathbf{u}_{2}-\mathbf{u}_{1})\cdot(\delta\mathbf{C}_{1}+\delta
C_{2}\times\mathbf{r})dx$
$\displaystyle+\rho_{f}\int_{B}(\mathbf{u}_{1}-\mathbf{u}_{2})\cdot\delta\mathbf{u}_{1}dx+2\mu\Delta
t\int_{B}\mathbf{D}(\mathbf{u}_{1}-\mathbf{u}_{2}):\mathbf{D}(\delta\mathbf{u}_{1})dx$
$\displaystyle=\rho_{f}\int_{B}\delta\mathbf{y}\cdot\boldsymbol{\pi}_{2}dx+\frac{\mu}{l_{s}}\int_{\gamma}\left(\Delta
t\delta\mathbf{u}_{1}-(\delta\mathbf{C}_{1}+\delta
C_{2}\times\mathbf{r})\right)\cdot\boldsymbol{\pi}_{2}d\gamma-\Delta
t\int_{B}\nabla\delta p_{2}\cdot\boldsymbol{\pi}_{2}dx$ $\displaystyle+\Delta
t\int_{B}\nabla
p_{\boldsymbol{\pi}_{2}}\cdot\left(\delta\mathbf{u}_{2}-\frac{1}{\Delta
t}(\delta\mathbf{C}_{1}+\delta
C_{2}\times\mathbf{r})\right)dx+\frac{\rho_{f}}{\Delta
t}\int_{B}(\mathbf{u}_{2}-\mathbf{u}_{1}-\boldsymbol{\pi}_{2})\cdot(\delta\mathbf{C}_{1}+\delta
C_{2}\times\mathbf{r})dx$
$\displaystyle+\rho_{f}\int_{B}(\mathbf{u}_{1}-\mathbf{u}_{2})\cdot\delta\mathbf{u}_{1}dx+2\mu\Delta
t\int_{B}\mathbf{D}(\mathbf{u}_{1}-\mathbf{u}_{2}):\mathbf{D}(\delta\mathbf{u}_{1})dx.$
(3.27)

If $q=\delta p_{2}$ in (3.26) and $q=p_{\boldsymbol{\pi}_{2}}$ in (3.24), we
have

$\displaystyle\delta
J(\mathbf{y},\mathbf{C}_{1},C_{2})=\rho_{f}\int_{B}\delta\mathbf{y}\cdot\boldsymbol{\pi}_{2}dx+\frac{\mu}{l_{s}}\int_{\gamma}\left(\Delta
t\delta\mathbf{u}_{1}-(\delta\mathbf{C}_{1}+\delta
C_{2}\times\mathbf{r})\right)\cdot\boldsymbol{\pi}_{2}d\gamma$
$\displaystyle+\frac{\rho_{f}}{\Delta
t}\int_{B}(\mathbf{u}_{2}-\mathbf{u}_{1}-\boldsymbol{\pi}_{2})\cdot\left(\delta\mathbf{C}_{1}+\delta
C_{2}\times\mathbf{r})\right)dx$
$\displaystyle+\rho_{f}\int_{B}(\mathbf{u}_{1}-\mathbf{u}_{2})\cdot\delta\mathbf{u}_{1}dx+2\mu\Delta
t\int_{B}\mathbf{D}(\mathbf{u}_{1}-\mathbf{u}_{2}):\mathbf{D}(\delta\mathbf{u}_{1})dx.$
(3.28)

We introduce now $\boldsymbol{\pi}_{1}$ and $p_{\boldsymbol{\pi}_{1}}$
solution of the following linear variational problem:

$\displaystyle\boldsymbol{\pi}_{1}\in(H^{1}_{0}(\Omega))^{2},\
p_{\boldsymbol{\pi}_{1}}\in H^{1}(\Omega),$
$\displaystyle\rho_{f}\int_{\Omega}\boldsymbol{\pi}_{1}\cdot\mathbf{v}dx+2\mu\Delta
t\int_{\Omega}\mathbf{D}(\boldsymbol{\pi}_{1}):\mathbf{D}(\mathbf{v})dx+\Delta
t\int_{\Omega}\nabla p_{\boldsymbol{\pi}_{1}}\cdot\mathbf{v}dx$ (3.29)
$\displaystyle=\rho_{f}\int_{B}(\mathbf{u}_{1}-\mathbf{u}_{2})\cdot\mathbf{v}dx+2\mu\Delta
t\int_{B}\mathbf{D}(\mathbf{u}_{1}-\mathbf{u}_{2})\cdot\mathbf{D}(\mathbf{v})dx+\frac{\mu\Delta
t}{l_{s}}\int_{\gamma}\boldsymbol{\pi}_{2}\cdot\mathbf{v}d\gamma,$
$\displaystyle\forall\mathbf{v}\in(H^{1}_{0}(\Omega))^{2},$
$\displaystyle\int_{\Omega}\boldsymbol{\pi}_{1}\cdot\nabla qdx=0,\ \forall
q\in H^{1}(\Omega).$ (3.30)

Taking $\mathbf{v}=\delta\mathbf{u}_{1}$ in (3.29) and combining with (3.28),
we obtain

$\displaystyle\delta
J(\mathbf{y},\mathbf{C}_{1},C_{2})=\rho_{f}\int_{B}\delta\mathbf{y}\cdot\boldsymbol{\pi}_{2}dx+\rho_{f}\int_{\Omega}\boldsymbol{\pi}_{1}\cdot\delta\mathbf{u}_{1}dx+2\mu\Delta
t\int_{\Omega}\mathbf{D}(\boldsymbol{\pi}_{1}):\mathbf{D}(\delta\mathbf{u}_{1})dx$
$\displaystyle+\Delta t\int_{\Omega}\nabla
p_{\boldsymbol{\pi}_{1}}\cdot\delta\mathbf{u}_{1}dx+\frac{\rho_{f}}{\Delta
t}\int_{B}(\mathbf{u}_{2}-\mathbf{u}_{1}-\boldsymbol{\pi}_{2})\cdot(\delta\mathbf{C}_{1}+\delta
C_{2}\times\mathbf{r})dx-\frac{\mu}{l_{s}}\int_{\gamma}(\delta\mathbf{C}_{1}+\delta
C_{2}\times\mathbf{r})\cdot\boldsymbol{\pi}_{2}d\gamma.$ (3.31)

If $q=p_{\boldsymbol{\pi}_{1}}$ in (3.22), we have

$\displaystyle\delta
J(\mathbf{y},\mathbf{C}_{1},C_{2})=\rho_{f}\int_{B}\delta\mathbf{y}\cdot\boldsymbol{\pi}_{2}dx+\rho_{f}\int_{\Omega}\boldsymbol{\pi}_{1}\cdot\delta\mathbf{u}_{1}dx+2\mu\Delta
t\int_{\Omega}\mathbf{D}(\boldsymbol{\pi}_{1}):\mathbf{D}(\delta\mathbf{u}_{1})dx$
$\displaystyle+\frac{\rho_{f}}{\Delta
t}\int_{B}(\mathbf{u}_{2}-\mathbf{u}_{1}-\boldsymbol{\pi}_{2})\cdot(\delta\mathbf{C}_{1}+\delta
C_{2}\times\mathbf{r})dx-\frac{\mu}{l_{s}}\int_{\gamma}(\delta\mathbf{C}_{1}+\delta
C_{2}\times\mathbf{r})\cdot\boldsymbol{\pi}_{2}d\gamma.$ (3.32)

Taking $\mathbf{v}=\boldsymbol{\pi}_{1}$ in (3.21) and combining with (3.32),
we finally obtain

$\displaystyle\delta
J(\mathbf{y},\mathbf{C}_{1},C_{2})=\rho_{f}\int_{B}\delta\mathbf{y}\cdot\boldsymbol{\pi}_{2}dx+\rho_{f}\int_{B}\delta\mathbf{y}\cdot\boldsymbol{\pi}_{1}dx-\Delta
t\int_{\Omega}\nabla\delta p_{1}\cdot\boldsymbol{\pi}_{1}dx$
$\displaystyle+\frac{\rho_{f}}{\Delta
t}\int_{B}(\mathbf{u}_{2}-\mathbf{u}_{1}-\boldsymbol{\pi}_{2})\cdot(\delta\mathbf{C}_{1}+\delta
C_{2}\times\mathbf{r})dx-\frac{\mu}{l_{s}}\int_{\gamma}(\delta\mathbf{C}_{1}+\delta
C_{2}\times\mathbf{r})\cdot\boldsymbol{\pi}_{2}d\gamma.$ (3.33)

If $q=\delta p_{1}$ in (3.30), we have

$\displaystyle\delta
J(\mathbf{y},\mathbf{C}_{1},C_{2})=\rho_{f}\int_{B}(\boldsymbol{\pi}_{1}+\boldsymbol{\pi}_{2})\cdot\delta\mathbf{y}dx+\frac{\rho_{f}}{\Delta
t}\int_{B}(\mathbf{u}_{2}-\mathbf{u}_{1}-\boldsymbol{\pi}_{2})\cdot(\delta\mathbf{C}_{1}+\delta
C_{2}\times\mathbf{r})dx$
$\displaystyle-\frac{\mu}{l_{s}}\int_{\gamma}(\delta\mathbf{C}_{1}+\delta
C_{2}\times\mathbf{r})\cdot\boldsymbol{\pi}_{2}d\gamma.$ (3.34)

Therefore,

$\displaystyle\frac{\partial
J}{\partial\mathbf{y}}(\mathbf{y},\mathbf{C}_{1},C_{2})=\rho_{f}\left(\boldsymbol{\pi}_{1}|_{B}+\boldsymbol{\pi}_{2}\right),$
(3.35) $\displaystyle\frac{\partial
J}{\partial\mathbf{C}_{1}}(\mathbf{y},\mathbf{C}_{1},C_{2})=\frac{\rho_{f}}{\Delta
t}\int_{B}(\mathbf{u}_{2}-\mathbf{u}_{1}-\boldsymbol{\pi}_{2})dx-\frac{\mu}{l_{s}}\int_{\gamma}\boldsymbol{\pi}_{2}d\gamma,$
(3.36) $\displaystyle\frac{\partial J}{\partial
C_{2}}(\mathbf{y},\mathbf{C}_{1},C_{2})=\frac{\rho_{f}}{\Delta
t}\int_{B}\left(-(\mathbf{u}_{2}-\mathbf{u}_{1}-\boldsymbol{\pi}_{2})_{x}\mathbf{r}_{y}+(\mathbf{u}_{2}-\mathbf{u}_{1}-\boldsymbol{\pi}_{2})_{y}\mathbf{r}_{x}\right)dx$
$\displaystyle-\frac{\mu}{l_{s}}\int_{\gamma}\left(-(\boldsymbol{\pi}_{2})_{x}\mathbf{r}_{y}+(\boldsymbol{\pi}_{2})_{y}\mathbf{r}_{x}\right)d\gamma,$
(3.37)

where $(\cdot)_{x}$ and $(\cdot)_{y}$ are the x-component and y-component,
respectively.

####  3.1.3 A conjugate gradient algorithm to the solution of the
least–squares problem (3.18)–(3.19)

Taking into account the way $\frac{\partial
J}{\partial\mathbf{y}}(\mathbf{y},\mathbf{C}_{1},C_{2})$, $\ \frac{\partial
J}{\partial\mathbf{C}_{1}}(\mathbf{y},\mathbf{C}_{1},C_{2})$ and
$\frac{\partial J}{\partial C_{2}}(\mathbf{y},\mathbf{C}_{1},C_{2})$, we are
going to give a more practical formulation. The algorithm reads as

$\displaystyle\mathbf{y}^{0},\ \mathbf{C}^{0}_{1},\ C_{2}^{0}\ \mbox{are given
in}\ (L^{2}(B))^{2}\times\mathbb{R}^{2}\times\mathbb{R},$ (3.38)

solve the following saddle point problems

$\displaystyle\mathbf{u}^{0}_{1}\in(H^{1}(\Omega))^{2},\
\mathbf{u}^{0}_{1}=\mathbf{g}_{\Gamma}\ \mbox{on}\ \Gamma,\ p^{0}_{1}\in
H^{1}(\Omega),$
$\displaystyle\rho_{f}\int_{\Omega}\mathbf{u}^{0}_{1}\cdot\mathbf{v}dx+2\mu\Delta
t\int_{\Omega}\mathbf{D}(\mathbf{u}^{0}_{1}):\mathbf{D}(\mathbf{v})dx+\Delta
t\int_{\Omega}\nabla p^{0}_{1}\cdot\mathbf{v}dx$
$\displaystyle=\rho_{f}\int_{B}\mathbf{y}^{0}\cdot\mathbf{v}dx+\rho_{f}\Delta
t\int_{\Omega\backslash\bar{B}}\mathbf{g}\cdot\mathbf{v}dx+\rho_{f}\int_{\Omega}\mathbf{u}_{\ast}\cdot\mathbf{v}dx,\
\forall\mathbf{v}\in(H_{0}^{1}(\Omega))^{2},$ (3.39)
$\displaystyle\int_{\Omega}\mathbf{u}^{0}_{1}\cdot\nabla
qdx=\int_{\Gamma}(\mathbf{g}_{\Gamma}\cdot\mathbf{n})qd\Gamma,\ \forall q\in
H^{1}(\Omega),$ (3.40) $\displaystyle\mathbf{u}^{0}_{2}\in(H^{1}(B))^{2},\
\left(\mathbf{u}^{0}_{2}-\left(\widetilde{\mathbf{U}}+\widetilde{\omega}\times\mathbf{r}\right)-\frac{1}{\Delta
t}\left(\mathbf{C}^{0}_{1}+C^{0}_{2}\times\mathbf{r}\right)\right)\cdot\mathbf{n}=0\
\mbox{on}\ \gamma,\ p^{0}_{2}\in H^{1}(B),$
$\displaystyle\rho_{f}\int_{B}\mathbf{u}^{0}_{2}\cdot\mathbf{v}dx+2\mu\Delta
t\int_{B}\mathbf{D}(\mathbf{u}^{0}_{2}):\mathbf{D}(\mathbf{v})dx+\Delta
t\int_{B}\nabla p^{0}_{2}\cdot\mathbf{v}dx$
$\displaystyle=\rho_{f}\int_{B}\mathbf{y}^{0}\cdot\mathbf{v}dx+\rho_{f}\int_{B}\mathbf{u}_{\ast}\cdot\mathbf{v}dx$
$\displaystyle+\frac{\mu}{l_{s}}\int_{\gamma}\left(\Delta
t\mathbf{u}^{0}_{1}-\Delta
t\left(\widetilde{\mathbf{U}}+\widetilde{\omega}\times\mathbf{r}\right)-\left(\mathbf{C}^{0}_{1}+C^{0}_{2}\times\mathbf{r}\right)\right)\cdot\mathbf{v}d\gamma,\
\forall\mathbf{v}\in\mathbb{V}_{0B},$ (3.41)
$\displaystyle\int_{B}\mathbf{u}^{0}_{2}\cdot\nabla
qdx=\int_{\gamma}\left(\left(\widetilde{\mathbf{U}}+\widetilde{\omega}\times\mathbf{r}\right)+\frac{1}{\Delta
t}\left(\mathbf{C}^{0}_{1}+C^{0}_{2}\times\mathbf{r}\right)\right)\cdot\mathbf{n}qd\gamma,\
\forall q\in H^{1}(B),$ (3.42)
$\displaystyle\boldsymbol{\pi}^{0}_{2}\in\mathbb{V}_{0B},\
p^{0}_{\boldsymbol{\pi}_{2}}\in H^{1}(B),$
$\displaystyle\rho_{f}\int_{B}\boldsymbol{\pi}^{0}_{2}\cdot\mathbf{v}dx+2\mu\Delta
t\int_{B}\mathbf{D}(\boldsymbol{\pi}^{0}_{2}):\mathbf{D}(\mathbf{v})dx+\Delta
t\int_{B}\nabla p^{0}_{\boldsymbol{\pi}_{2}}\cdot\mathbf{v}dx$ (3.43)
$\displaystyle=\rho_{f}\int_{B}(\mathbf{u}^{0}_{2}-\mathbf{u}^{0}_{1})\cdot\mathbf{v}dx+2\mu\Delta
t\int_{B}\mathbf{D}(\mathbf{u}^{0}_{2}-\mathbf{u}^{0}_{1}):\mathbf{D}(\mathbf{v}),\forall\mathbf{v}\in\mathbb{V}_{0B},$
$\displaystyle\int_{B}\boldsymbol{\pi}^{0}_{2}\cdot\nabla qdx=0,\ \forall q\in
H^{1}(B),$ (3.44)
$\displaystyle\boldsymbol{\pi}^{0}_{1}\in(H^{1}_{0}(\Omega))^{2},\
p^{0}_{\boldsymbol{\pi}_{1}}\in H^{1}(\Omega),$
$\displaystyle\rho_{f}\int_{\Omega}\boldsymbol{\pi}^{0}_{1}\cdot\mathbf{v}dx+2\mu\Delta
t\int_{\Omega}\mathbf{D}(\boldsymbol{\pi}^{0}_{1}):\mathbf{D}(\mathbf{v})dx+\Delta
t\int_{\Omega}\nabla p^{0}_{\boldsymbol{\pi}_{1}}\cdot\mathbf{v}dx$
$\displaystyle=\rho_{f}\int_{B}(\mathbf{u}^{0}_{1}-\mathbf{u}^{0}_{2})\cdot\mathbf{v}dx+2\mu\Delta
t\int_{B}\mathbf{D}(\mathbf{u}^{0}_{1}-\mathbf{u}^{0}_{2}):\mathbf{D}(\mathbf{v})+\frac{\mu}{l_{s}}\Delta
t\int_{\gamma}\boldsymbol{\pi}^{0}_{2}\cdot\mathbf{v}d\gamma,$
$\displaystyle\forall\mathbf{v}\in(H^{1}_{0}(\Omega))^{2},$ (3.45)
$\displaystyle\int_{\Omega}\boldsymbol{\pi}^{0}_{1}\cdot\nabla qdx=0,\ \forall
q\in H^{1}(\Omega),$ (3.46)

and set

$\displaystyle\mathbf{g}^{0}_{1}=\rho_{f}(\boldsymbol{\pi}^{0}_{1}|_{B}+\boldsymbol{\pi}^{0}_{2}),$
$\displaystyle\mathbf{w}^{0}_{1}=\mathbf{g}^{0}_{1},$
$\displaystyle\mathbf{g}^{0}_{2}=\frac{\rho_{f}}{\Delta
t}\int_{B}(\mathbf{u}^{0}_{2}-\mathbf{u}^{0}_{1}-\boldsymbol{\pi}^{0}_{2})dx-\frac{\mu}{l_{s}}\int_{\gamma}\boldsymbol{\pi}^{0}_{2}d\gamma,$
$\displaystyle\mathbf{w}^{0}_{2}=\mathbf{g}^{0}_{2},$ $\displaystyle
g^{0}_{3}=\frac{\rho_{f}}{\Delta
t}\int_{B}\left(-(\mathbf{u}^{0}_{2}-\mathbf{u}^{0}_{1}-\boldsymbol{\pi}^{0}_{2})_{x}\mathbf{r}_{y}+(\mathbf{u}^{0}_{2}-\mathbf{u}^{0}_{1}-\boldsymbol{\pi}^{0}_{2})_{y}\mathbf{r}_{x}\right)dx$
$\displaystyle-\frac{\mu}{l_{s}}\int_{\gamma}\left(-(\boldsymbol{\pi}^{0}_{2})_{x}\mathbf{r}_{y}+(\boldsymbol{\pi}^{0}_{2})_{y}\mathbf{r}_{x}\right)d\gamma,$
$\displaystyle w^{0}_{3}=g^{0}_{3}.$

For $k\geq
0,\mathbf{y}^{k},\mathbf{C}^{k}_{1},C^{k}_{2},\mathbf{g}^{k}_{1},\mathbf{g}^{k}_{2},g^{k}_{3}$
and $\mathbf{w}^{k}_{1},\mathbf{w}^{k}_{2},w^{k}_{3}$ being known, we compute
$\mathbf{y}^{k+1},\mathbf{C}^{k+1}_{1},C^{k+1}_{2}$,
$\mathbf{g}^{k+1}_{1},\mathbf{g}^{k+1}_{2},g_{3}^{k+1}$ and
$\mathbf{w}_{1}^{k+1},\mathbf{w}_{2}^{k+1},w_{3}^{k+1}$ as follows:
solve the saddle point problems

$\displaystyle\overline{\mathbf{u}}^{k}_{1}\in(H^{1}_{0}(\Omega))^{2},\
\overline{p}^{k}_{1}\in H^{1}(\Omega),$
$\displaystyle\rho_{f}\int_{\Omega}\overline{\mathbf{u}}^{k}_{1}\cdot\mathbf{v}dx+2\mu\Delta
t\int_{\Omega}\mathbf{D}(\overline{\mathbf{u}}^{k}_{1}):\mathbf{D}(\mathbf{v})dx+\Delta
t\int_{\Omega}\nabla\overline{p}^{k}_{1}\cdot\mathbf{v}dx$
$\displaystyle=\rho_{f}\int_{B}\mathbf{w}^{k}_{1}\cdot\mathbf{v}dx,\
\forall\mathbf{v}\in(H_{0}^{1}(\Omega))^{2},$ (3.47)
$\displaystyle\int_{\Omega}\overline{\mathbf{u}}^{k}_{1}\cdot\nabla qdx=0,\
\forall q\in H^{1}(\Omega),$ (3.48)
$\displaystyle\overline{\mathbf{u}}^{k}_{2}\in\mathbb{V}_{0B},\
\overline{p}^{k}_{2}\in H^{1}(B),$
$\displaystyle\rho_{f}\int_{B}\overline{\mathbf{u}}^{k}_{2}\cdot\mathbf{v}dx+2\mu\Delta
t\int_{B}\mathbf{D}(\overline{\mathbf{u}}^{k}_{2}):\mathbf{D}(\mathbf{v})dx+\Delta
t\int_{B}\nabla\overline{p}^{k}_{2}\cdot\mathbf{v}dx$
$\displaystyle=\rho_{f}\int_{B}\mathbf{w}^{k}\cdot\mathbf{v}dx+\frac{\mu}{l_{s}}\int_{\gamma}\left(\Delta
t\overline{\mathbf{u}}^{k}_{1}-\left(\mathbf{w}^{k}_{2}+w^{k}_{3}\times\mathbf{r}\right)\right)\cdot\mathbf{v}d\gamma,\
\forall\mathbf{v}\in\mathbb{V}_{0B},$ (3.49)
$\displaystyle\int_{B}\overline{\mathbf{u}}^{k}_{2}\cdot\nabla qdx=0,\ \forall
q\in H^{1}(B),$ (3.50)
$\displaystyle\overline{\boldsymbol{\pi}}^{k}_{2}\in\mathbb{V}_{0B},\
\overline{p}^{k}_{\boldsymbol{\pi}_{2}}\in H^{1}(B),$
$\displaystyle\rho_{f}\int_{B}\overline{\boldsymbol{\pi}}^{k}_{2}\cdot\mathbf{v}dx+2\mu\Delta
t\int_{B}\mathbf{D}(\overline{\boldsymbol{\pi}}^{k}_{2}):\mathbf{D}(\mathbf{v})dx+\Delta
t\int_{B}\nabla\overline{p}^{k}_{\boldsymbol{\pi}_{2}}\cdot\mathbf{v}dx$
(3.51)
$\displaystyle=\rho_{f}\int_{B}(\overline{\mathbf{u}}^{k}_{2}-\overline{\mathbf{u}}^{k}_{1})\cdot\mathbf{v}dx+2\mu\Delta
t\int_{B}\mathbf{D}(\overline{\mathbf{u}}^{k}_{2}-\overline{\mathbf{u}}^{k}_{1}):\mathbf{D}(\mathbf{v}),\forall\mathbf{v}\in\mathbb{V}_{0B},$
$\displaystyle\int_{B}\overline{\boldsymbol{\pi}}^{k}_{2}\cdot\nabla
qdx=0,\forall q\in H^{1}(B),$ (3.52)
$\displaystyle\overline{\boldsymbol{\pi}}^{k}_{1}\in(H^{1}_{0}(\Omega))^{2},\
\overline{p}^{k}_{\boldsymbol{\pi}_{1}}\in H^{1}(\Omega),$
$\displaystyle\rho_{f}\int_{\Omega}\overline{\boldsymbol{\pi}}^{k}_{1}\cdot\mathbf{v}dx+2\mu\Delta
t\int_{\Omega}\mathbf{D}(\overline{\boldsymbol{\pi}}^{k}_{1}):\mathbf{D}(\mathbf{v})dx+\Delta
t\int_{\Omega}\nabla\overline{p}^{k}_{\boldsymbol{\pi}_{1}}\cdot\mathbf{v}dx$
$\displaystyle=\rho_{f}\int_{B}(\overline{\mathbf{u}}^{k}_{1}-\overline{\mathbf{u}}^{k}_{2})\cdot\mathbf{v}dx+2\mu\Delta
t\int_{B}\mathbf{D}(\overline{\mathbf{u}}^{k}_{1}-\overline{\mathbf{u}}^{k}_{2}):\mathbf{D}(\mathbf{v})+\frac{\mu}{l_{s}}\Delta
t\int_{\gamma}\overline{\boldsymbol{\pi}}^{k}_{2}\cdot\mathbf{v}d\gamma,$
$\displaystyle\forall\mathbf{v}\in(H^{1}_{0}(\Omega))^{2},$ (3.53)
$\displaystyle\int_{\Omega}\overline{\boldsymbol{\pi}}^{k}_{1}\cdot\nabla
qdx=0,\ \forall q\in H^{1}(\Omega),$ (3.54)

and set

$\displaystyle\overline{\mathbf{g}}^{k}_{1}=\rho_{f}(\overline{\boldsymbol{\pi}}^{k}_{1}|_{B}+\bar{\boldsymbol{\pi}}^{k}_{2}),$
(3.55) $\displaystyle\overline{\mathbf{g}}^{k}_{2}=\frac{\rho_{f}}{\Delta
t}\int_{B}(\overline{\mathbf{u}}^{k}_{2}-\overline{\mathbf{u}}^{k}_{1}-\overline{\boldsymbol{\pi}}^{k}_{2})dx-\frac{\mu}{l_{s}}\int_{\gamma}\overline{\boldsymbol{\pi}}^{k}_{2}d\gamma,$
(3.56) $\displaystyle\overline{g}^{k}_{3}=\frac{\rho_{f}}{\Delta
t}\int_{B}\left(-\left(\overline{\mathbf{u}}^{k}_{2}-\overline{\mathbf{u}}^{k}_{1}-\overline{\boldsymbol{\pi}}^{k}_{2}\right)_{x}\mathbf{r}_{y}+(\overline{\mathbf{u}}^{k}_{2}-\overline{\mathbf{u}}^{k}_{1}-\overline{\boldsymbol{\pi}}^{k}_{2})_{y}\mathbf{r}_{x}\right)dx$
$\displaystyle-\frac{\mu}{l_{s}}\int_{\gamma}\left(-(\overline{\boldsymbol{\pi}}^{k}_{2})_{x}\mathbf{r}_{y}+(\overline{\boldsymbol{\pi}}^{k}_{2})_{y}\mathbf{r}_{x}\right)d\gamma.$
(3.57)

Compute

$\displaystyle\lambda_{k}=\frac{\int_{B}|\mathbf{g}^{k}_{1}|^{2}dx+|\mathbf{g}^{k}_{2}|^{2}+|g^{k}_{3}|^{2}}{\int_{B}\overline{\mathbf{g}}^{k}_{1}\cdot\mathbf{w}^{k}_{1}dx+\overline{\mathbf{g}}^{k}_{2}\cdot\mathbf{w}^{k}_{2}+\overline{g}^{k}_{3}w^{k}_{3}},$
(3.58)
$\displaystyle\mathbf{y}^{k+1}=\mathbf{y}^{k}-\lambda_{k}\mathbf{w}^{k}_{1},$
(3.59)
$\displaystyle\mathbf{C}^{k+1}_{1}=\mathbf{C}^{k}_{1}-\lambda_{k}\mathbf{w}^{k}_{2},$
(3.60) $\displaystyle C^{k+1}_{2}=C^{k}_{2}-\lambda_{k}w^{k}_{3},$ (3.61)
$\displaystyle\mathbf{g}^{k+1}_{1}=\mathbf{g}^{k}_{1}-\lambda_{k}\overline{\mathbf{g}}^{k}_{1},$
(3.62)
$\displaystyle\mathbf{g}^{k+1}_{2}=\mathbf{g}^{k}_{2}-\lambda_{k}\overline{\mathbf{g}}^{k}_{2},$
(3.63) $\displaystyle g^{k+1}_{3}=g^{k}_{3}-\lambda_{k}\overline{g}^{k}_{3}.$
(3.64)

If
$\frac{\int_{B}|\mathbf{g}^{k+1}_{1}|^{2}dx+|\mathbf{g}^{k+1}_{2}|^{2}+|g^{k+1}_{3}|^{2}}{\int_{B}|\mathbf{g}^{0}_{1}|^{2}dx+|\mathbf{g}^{0}_{2}|^{2}+|g^{0}_{3}|^{2}}\leq\
\mbox{tol}$, we take $\mathbf{y}=\mathbf{y}^{k+1},\
\mathbf{C}_{1}=\mathbf{C}^{k+1}_{1},$ and $C_{2}=C_{2}^{k+1}$, which is
substituted into (3.10)–(3.11) and (3.14)–(3.15), we get

$\displaystyle\mathbf{u}=\mathbf{u}_{1}|_{\Omega\backslash\overline{B}},$
$\displaystyle\mathbf{U}=\widetilde{\mathbf{U}}+\frac{1}{\Delta
t}\mathbf{C}_{1},$ $\displaystyle\omega=\widetilde{\omega}+\frac{1}{\Delta
t}C_{2}.$

Otherwise, compute

$\displaystyle\gamma_{k}=\frac{\int_{B}|\mathbf{g}^{k+1}_{1}|^{2}dx+|\mathbf{g}^{k+1}_{2}|^{2}+|g^{k+1}_{3}|^{2}}{\int_{B}|\mathbf{g}^{k}_{1}|^{2}dx+|\mathbf{g}^{k}_{2}|^{2}+|g^{k}_{3}|^{2}},$
(3.65)

and set

$\displaystyle\mathbf{w}^{k+1}_{1}=\mathbf{g}^{k+1}_{1}+\gamma_{k}\mathbf{w}^{k}_{1},$
(3.66)
$\displaystyle\mathbf{w}^{k+1}_{2}=\mathbf{g}^{k+1}_{2}+\gamma_{k}\mathbf{w}^{k}_{2},$
(3.67) $\displaystyle w^{k+1}_{3}=g^{k+1}_{3}+\gamma_{k}w^{k}_{3}.$ (3.68)

Do $n+1\rightarrow n$ and return to (3.47).

By inspection of (3.39)–(3.40), (3.45)–(3.46), (3.47)–(3.48) and
(3.53)–(3.54), it is clear that the pressure term can be treated as a
corresponding Lagrange multiplier in the space $H_{0}$, which is defined as
$H_{0}=\\{q|q\in H^{1}(\Omega),\ \int_{\Omega}qdx=0\\}$. Similar arguments can
be applied to (3.41)–(3.42), (3.43)–(3.44), (3.49)–(3.50) and (3.51)–(3.52).
Therefore, a preconditioned conjugate gradient algorithm is easy to implement
and seems to have excellent convergence properties.

######  Remark 3.2.

If $\mathbf{y}^{0}$, $\mathbf{C}^{0}_{1}$, $C^{0}_{2}$ are close to
$\mathbf{y}$, $\mathbf{C}_{1}$, $C_{2}$ (the stopping criterion we used for
algorithm (3.38)–(3.68)) may lead to more iterations than necessary. A more
realistic stopping test is given by

$\displaystyle\frac{\int_{B}|\mathbf{g}^{n+1}_{1}|^{2}dx+|\mathbf{g}^{k+1}_{2}|^{2}+|g^{k+1}_{3}|^{2}}{\max\left(\int_{B}|\mathbf{g}^{0}_{1}|^{2}dx+|\mathbf{g}^{0}_{2}|^{2}+|g^{0}_{3}|^{2},\int_{B}|\mathbf{y}^{n+1}|^{2}dx+|\mathbf{C}^{k+1}_{1}|^{2}+|C^{k+1}_{2}|^{2}\right)}\leq\
\mbox{tol}.$

Other stopping criteria can be used.

##  4 Finite element realization of the operator-splitting scheme
(2.14)–(2.18) and (2.8)–(2.13)

###  4.1 Generalities

We describe in this section the finite element implementation of the
least–squares/ fictitious domain methodology discussed in above sections. We
will assume that $\overline{B}\subset\Omega\subset\mathbb{R}^{2}$ and that
$\Omega$ is convex and/or has a smooth boundary; similarly, we assume that
$\gamma$ is smooth. Our approximation of choice will be the
_Bercovier–Pironneau_ one, discussed with details in Glowinski[2](see also the
references therein). For simplicity we still denote by $\Omega$ and $B$ the
polygonal approximations of the above domains. From the triangulations
$\mathcal{T}_{h_{1}}$ of $\Omega$ and $\mathcal{T}_{h_{2}}$ of $B$ we define
the following finite dimensional spaces. The pressure space $P_{h_{1}}$ and
$P_{Bh_{2}}$ are

$\displaystyle P_{h_{1}}=\\{q|q\in C^{0}(\overline{\Omega}),q|_{T}\in
P_{1},\forall T\in\mathcal{T}_{h_{1}}\\},$ (4.1) $\displaystyle
P_{Bh_{2}}=\\{q|q\in C^{0}(\overline{B}),q|_{T}\in P_{1},\forall
T\in\mathcal{T}_{h_{2}}\\},$ (4.2) $\displaystyle P_{0h_{1}}=\\{q|q\in
P_{h_{1}},\ \int_{\Omega}qdx=0\\},$ (4.3) $\displaystyle P_{0Bh_{2}}=\\{q|q\in
P_{Bh_{2}},\ \int_{B}qdx=0\\},$ (4.4)

where $P_{1}$ being the space of the polynomials of two variables of degree
$\leq 1$ and $h_{1}$ (resp., $h_{2}$) the length of the largest edge(s) of the
finite element triangulation $\mathcal{T}_{h_{1}}$ (resp.,
$\mathcal{T}_{h_{2}}$). Next, we divide each triangle of
$\mathcal{T}_{h_{1}}$(resp. $\mathcal{T}_{h_{2}}$) in four sub-triangles, by
joining the mid-points and denote by $\mathcal{T}_{h_{1}/2}$(resp.
$\mathcal{T}_{h_{2}/2}$) the resulting triangulation. From
$\mathcal{T}_{h_{1}/2}$(resp. $\mathcal{T}_{h_{2}/2}$), we define the velocity
spaces

$\displaystyle\mathbb{V}_{h_{1}}=\\{\mathbf{v}|\mathbf{v}\in
C^{0}(\overline{\Omega}),\mathbf{v}|_{T}\in(P_{1})^{2},\forall
T\in\mathcal{T}_{h_{1}/2}\\},$ (4.5)
$\displaystyle\mathbb{V}_{0h_{1}}=\\{\mathbf{v}|\mathbf{v}\in\mathbb{V}_{h_{1}},\mathbf{v}=\mathbf{0}\
\mbox{on}\ \Gamma\\},$ (4.6)

and

$\displaystyle\mathbb{V}_{Bh_{2}}=\\{\mathbf{v}|\mathbf{v}\in
C^{0}(\overline{B}),\mathbf{v}|_{T}\in(P_{1})^{2},\forall
T\in\mathcal{T}_{h_{2}/2}\\}.$ (4.7)

The finite dimensional spaces $\mathbb{V}_{h_{1}},\mathbb{V}_{0h_{1}}$ and
$\mathbb{V}_{Bh_{2}}$ are finite dimensional approximations to
$(H^{1}(\Omega))^{2},(H_{0}^{1}(\Omega))^{2}$ and $(H^{1}(B))^{2}$,
respectively. Similarly, we will use $\mathbb{V}_{Bh_{2}}$ to approximate the
control space $(L^{2}(B))^{2}$. We use $\Sigma_{h_{1}}$ to denote the set of
the vertices of $\mathcal{T}_{h_{1}/2}$ and $\Sigma_{h_{2}}$ to denote the set
of the vertices of $\mathcal{T}_{h_{2}/2}$.

###  4.2 Finite element approximation of the least–squares problem
(3.18)–(3.19)

To approximate the least–squares problem (3.18)–(3.19), we suggest

$\displaystyle\mathbf{y}^{\ast}_{h}\in\mathbb{V}_{Bh_{2}},\
\mathbf{C}^{\ast}_{1}\in\mathbb{R}^{2},\ C^{\ast}_{2}\in\mathbb{R},$
$\displaystyle
J_{h}(\mathbf{y}^{\ast}_{h},\mathbf{C}^{\ast}_{1},C^{\ast}_{2})\leq
J_{h}(\mathbf{y}_{h},\mathbf{C}_{1},C_{2}),\
\forall\mathbf{y}_{h}\in\mathbb{V}_{Bh_{2}},\
\mathbf{C}_{1}\in\mathbb{R}^{2},\ C_{2}\in\mathbb{R}$ (4.8)

where

$\displaystyle
J_{h}(\mathbf{y}_{h},\mathbf{C}_{1},C_{2})=\frac{1}{2}\left[\rho_{f}\int_{B}|\mathbf{u}_{2}-\Pi_{2}\mathbf{u}_{1}|^{2}dx+2\mu\Delta
t\int_{B}|\mathbf{D}(\mathbf{u}_{2}-\Pi_{2}\mathbf{u}_{1})|^{2}dx\right].$
(4.9)

In the above, $\mathbf{u}_{1}$ is solution to the following fully discrete
approximate saddle point problem:

$\displaystyle\mathbf{u}_{1}\in\mathbb{V}_{h_{1}},\
\mathbf{u}_{1}=\mathbf{g}_{\Gamma}\ \mbox{on}\ \Gamma,\ p_{1}\in P_{h_{1}},$
$\displaystyle\rho_{f}\int_{\Omega}\mathbf{u}_{1}\cdot\mathbf{v}dx+2\mu\Delta
t\int_{\Omega}\mathbf{D}(\mathbf{u}_{1}):\mathbf{D}(\mathbf{v})dx+\Delta
t\int_{\Omega}\nabla p_{1}\cdot\mathbf{v}dx$
$\displaystyle=\rho_{f}\int_{B}\mathbf{y}_{h}\cdot\Pi_{2}\mathbf{v}dx+\rho_{f}\int_{\Omega}\mathbf{u}_{\ast}\cdot\mathbf{v}dx+\rho_{f}\Delta
t\int_{\Omega\backslash\bar{B}}\mathbf{g}\cdot\mathbf{v}dx,\
\forall\mathbf{v}\in\mathbb{V}_{0h_{1}},$ (4.10)
$\displaystyle\int_{\Omega}\mathbf{u}_{1}\cdot\nabla
qdx=\int_{\Gamma}(\mathbf{g}_{\Gamma}\cdot\mathbf{n})qd\Gamma,\ \forall q\in
P_{h_{1}},$ (4.11)

where $\Pi_{2}:(C^{0}(\overline{\Omega}))^{2}\rightarrow\mathbb{V}_{h_{2}}$ is
the interpolation operator defined as follows

$\displaystyle\Pi_{2}\mathbf{v}=\sum_{i=1}^{N_{h_{2}}}\mathbf{v}(Y_{i})w_{2i},\
\forall\ \mathbf{v}\in(C^{0}(\overline{\Omega}))^{2},$ (4.12)

where $\\{Y_{i}\\}_{i=1}^{N_{h_{2}}}$ being the set of the vertices of
$\mathcal{T}_{h_{2}/2}$ and $w_{2i}$ the $P_{1}$\- shape function in
$\mathbb{V}_{h_{2}}$ associated with the vertex $Y_{i}$. Return to (4.9), the
function $\mathbf{u}_{2}$ is the solution of the following saddle point
problem

$\displaystyle\mathbf{u}_{2}\in\mathbb{V}_{Bh_{2}},\
\left(\mathbf{u}_{2}-(\widetilde{\mathbf{U}}+\widetilde{\omega}\times\mathbf{r})-\frac{1}{\Delta
t}\left(\mathbf{C}_{1}+C_{2}\times\mathbf{r}\right)\right)\cdot\mathbf{n}=0\
\mbox{on}\ \gamma,\ p_{2}\in P_{B_{h_{2}}},$
$\displaystyle\rho_{f}\int_{B}\mathbf{u}_{2}\cdot\mathbf{v}dx+2\mu\Delta
t\int_{B}\mathbf{D}(\mathbf{u}_{2}):\mathbf{D}(\mathbf{v})dx+\Delta
t\int_{B}\nabla p_{2}\cdot\mathbf{v}dx$ (4.13)
$\displaystyle=\rho_{f}\int_{B}\mathbf{y}_{h}\cdot\mathbf{v}dx+\frac{\mu\Delta
t}{l_{s}}\int_{\gamma}(\Pi_{2}\mathbf{u}_{1}-(\widetilde{\mathbf{U}}+\widetilde{\omega}\times\mathbf{r}))\cdot\mathbf{v}d\gamma-\frac{\mu}{l_{s}}\int_{\gamma}(\mathbf{C}_{1}+C_{2}\times\mathbf{r})\cdot\mathbf{v}d\gamma$
$\displaystyle+\rho_{f}\int_{B}\Pi_{2}\mathbf{u}_{\ast}\cdot\mathbf{v}dx,\forall\mathbf{v}\in\mathbb{V}_{0Bh_{2}},$
$\displaystyle\int_{B}\mathbf{u}_{2}\cdot\nabla
qdx=\int_{\gamma}\left(\left(\widetilde{\mathbf{U}}+\widetilde{\omega}\times\mathbf{r}\right)+\frac{1}{\Delta
t}\left(\mathbf{C}_{1}+C_{2}\times\mathbf{r}\right)\right)\cdot qd\gamma\
\forall q\in P_{Bh_{2}},$ (4.14)
$\displaystyle\mathbf{U}=\widetilde{\mathbf{U}}+\frac{1}{\Delta
t}\mathbf{C}_{1},$ (4.15)
$\displaystyle\omega=\widetilde{\omega}+\frac{1}{\Delta t}C_{2},$ (4.16)

where
$\mathbb{V}_{0Bh_{2}}\triangleq\\{\mathbf{v}|\mathbf{v}\in\mathbb{V}_{Bh_{2}},\
\mathbf{v}\cdot\mathbf{n}=0\ \mbox{on}\
\gamma\\},\widetilde{\mathbf{U}}=\mathbf{U}_{\ast}-\frac{\Delta
t}{\mathbf{M}}\int_{\gamma}\left(-\Pi_{2}p_{1\ast}\mathbf{E}+2\mu\mathbf{D}(\Pi_{2}\mathbf{u}_{\ast})\right)d\gamma,$
and $\widetilde{\omega}=\omega_{\ast}-\frac{\Delta
t}{\mathbf{I}}\int_{\gamma}\mathbf{r}\times\left(-\Pi_{2}p_{1\ast}\mathbf{E}+2\mu\mathbf{D}(\Pi_{2}\mathbf{u}_{\ast})\right)d\gamma$.

With the introduce of interpolation operator $\Pi_{2}$, a finite element
discretization of the conjugate gradient algorithms (3.38)–(3.68) is easy to
realize. As mentioned above, the finite element discretization of
(3.39)–(3.40), (3.45)–(3.46) , (3.47)–(3.48) and (3.53)–(3.54) can be easily
solved by a preconditioned conjugate gradient algorithm which is detailed
discussed in Glowinski[2]. The finite element discretization of (3.41)–(3.42),
(3.43)–(3.44), (3.49)–(3.50) and (3.51)–(3.52) also have been discribed in
[10].

##  5 Numerical experiments

###  5.1 Jeffery orbit

In Jeffery [12], an elliptic particle was positioned in a symmetric domain
filled with simple shear rate $\dot{\gamma}$. For Stokes flows, i.e. in the
absence of inertia, assuming there is no-slip between particle and fluid, the
dynamic equation of Jeffery orbit is given by

$\displaystyle\omega=\frac{\dot{\gamma}}{2}\left(-1+p\cos 2\Theta\right),$
(5.1)

therefore the elliptic particle undergoes a periodic tumbling, where $\omega$
is the angular velocity of the ellipse, $\Theta$ is the angular orientation,
$p$ is a parameter measuring the anisotropy of the ellipse, i.e.
$p=\frac{1-e^{2}}{1+e^{2}}$, where $e$ is the aspect ratio defined by
$e=\frac{b}{a}\leq 1$, $a$ and $b$ are the major radius and minor radius of
the elliptic particle respectively. Since the particle is positioned at the
center of a symmetric domain, it has zero translational velocity.

We now use the method described in preceding sections to numerically simulate
the interactions between fluid and rigid body. The ellipse rigid particle,
which is located at the domain center $(1,1)$, with $a=\frac{1}{4}$ and
$b=\frac{1}{8}$. The computational domain is $[0,2]\times[0,2]$. The solid
walls move with speed $W=2$, in $+x$ direction on the top and $-x$ direction
on the bottom, and the particle rotates with an angular velocity $\omega$
which is solved by the numerical scheme introduced in the above sections. The
theoretic shear rate $\dot{\gamma}=2$ and the anisotropy parameter $p=0.6$ can
be calculated from the given parameters. The external force is $\mathbf{g}=0$
and we choose $\rho_{f}=1$, the density of the ellipse particle is the same as
the fluid and $\mu=1.0$, instead of Navier–Stokes equations, we numerically
simulate Stokes equations. The first issue is to verify that $\mathbf{C}_{1}$
and $C_{2}$ will approach to $0$ when $\Delta t$ and $h_{1}$ converge to $0$,
which is due to that these two terms are correction terms for translation
velocity $\mathbf{U}$ and angular velocity $\omega$. The values of
$\mathbf{C}_{1}$ and $C_{2}$ for different $\Delta t$ are shown in Table 1,
where $\Delta t=O(h^{2}_{1})$. As suggested by [4] and [10], we take
$h_{2}=h_{1}$ in our numerical simulations.

We will address how angular velocity $\omega$ changes with the slip length
$l_{s}$. We will numerically verify the effective anisotropy of the particle
could be enhanced by the boundary slip, which is presented in [19]. Figure 2
shows snapshots of velocity field in different time, and the rotation is
observed and translation is almost zero. Next, we will verify the Jeffery
orbit formulation (5.1). Using our numerical schemes with slip length
$l_{s}=\frac{1}{20}$, we see the dependence of the angular velocity $\omega$
and the angular orientation $\Theta$ in Figure 3 with different mesh size,
which shows the convergence of our numerical method. The numerical errors
(oscillations) can be reduced by mesh refinement. In order to further
investigate the Jeffery obit theory with boundary slip, we plot the fitting
curve according to (5.1) and our numerical results in Figure 4, where the
fitted shear rate $\dot{\gamma}=1.97$ (close to theoretical value
$\dot{\gamma}=2$), the fitted anisotropy parameter $p=0.69$ , it is observed
that our numerical results show a good cosine curve. In Figure 5, we see that
as slip length increases, the anisotropy parameter $p$ is also increased,
which has been explained in [19] that the effective anisotropy of the particle
could be enhanced by the boundary slip. The numerical solution of the
anisotropy parameter $p$ is computed by equation (5.1), since numerical
solution of $\omega$ can be obtained from our numerical scheme. With different
slip length $l_{s}=0,\ 1/200,\ 1/100,\ 1/40,\ 1/20$, we plot the numerical
results of $\omega$ along with angular orientation $\Theta$ in Figure 6. It is
clear to see that our numerical method can resolve the effect of slip length
to the solution efficiently. It is proposed that the effective shape of the
particle has the major and minor axes given by $a-l_{s}$ and $b-l_{s}$ in
[19]. It is observed that the slip length affect the amplitude of
trigonometric function, which can be obtained by varying the value of axial
ratio $p$. The results are shown in Figure 6, which are fairly matched with
the that in [19].

Table 1:  The variations of $\mathbf{C}_{1}$ and $C_{2}$. | $\Delta t$ | $\mathbf{C}_{1}$ | $C_{2}$
---|---|---|---
| 2.0e-3 | 1.85e-3, 2.12e-3 | 2.68e-3
| 4.88e-4 | 4.26e-4, 6.35e-4 | 8.76e-4
| 1.22e-4 | 7.32e-5, 8.11e-5 | 1.05e-4
| 3.51e-5 | 8.34e-6, 9.23e-6 | 3.24e-5

Figure 2: Snapshots of velocity field when $t=0.01,\ 0.37,\ 0.77$ and $1.17$
for $l_{s}=1/20$. Figure 3: Comparison of numerical results of different mesh
size $h_{1}$ for a fixed $h_{2}$ with $l_{s}=1/20$. Figure 4: The angular
velocity $\omega$ is plotted as a function of the ellipse orientation $\Theta$
for $l_{s}=1/20$. Figure 5: Numerical solutions of $p$ versus to slip length
$l_{s}$. Figure 6: Numerical solutions of $\omega$ versus to $\Theta$ with
$l_{s}=0,\ 1/200,\ 1/100,\ 1/40,\ 1/20$.

###  5.2 Sedimentation of an elliptic particle in a viscous flow

We consider the sedimentation of an elliptic particle in a rectangular domain.
The ellipse particle is initially static and falls down driven by the gravity
$\mathbf{g}=(0,-9.8)^{T}$ and $\mu=0.1$. The $\rho_{f}=1.0$ and
$\rho_{s}=1.01$ are the density of fluid and solid, respectively. The
computational domain $\Omega=[0,0.5]\times[0,10]$. The elliptic particle is
initially located at coordinate $(0.25,9.7)$. The major radius $a$ and minor
radius $b$ of the elliptic particle are $0.0625$ and $0.03125$, respectively.
From the physical experiments and numerical simulations, we know that the
particle will undergo a motion called “drafting and tumbling”, which was first
numerically demonstrated in [1]. The fully incompresible Navier–Stokes are
considered. We run the simulations with different slip lengths
$ls=2d,d,0.4d,0.2d,0$, where $d=2a$. The instantaneous transverse coordinate
and longitudinal coordinate of the center of the ellipse are shown in Figure
7–8. The angular velocity is shown in Figure 9. From those results, we clearly
find that the particle with different boundary conditions initially perform
differential dynamic process and finally reach similarly pseudo-steady state.
It is observed that slip length will increase the drag coefficient, which can
be verified by the numerical results obtained by immersed boundary–lattice
Boltzmann method in [18]. Therefore, the particle falls more faster with
increased slip length.

Figure 7: Transverse coordinate of the ellipse center for $l_{s}=2d,\ 1d,\
0.4d,\ 0.2d,0$, $d=2a$.  Figure 8: Longitudinal coordinate of the ellipse
center for $l_{s}=2d,\ 1d,\ 0.4d,\ 0.2d,0$, $d=2a$. Figure 9: Angular velocity
the ellipse center for $l_{s}=2d,\ 1d,\ 0.4d,\ 0.2d,0$, $d=2a$.

##  6 Conclusion

In this article, we develop a least–squares/fictitious domain method for
directly numerical simulation of fluid particle interaction with Navier slip
boundary condition at the fluid–particle interface. The method relies on a
least–squares formulation of the virtual control type, making it well–suited
for solution by a conjugate gradient algorithm operating in a well–chosen
control space. Numerical experiments show that our method works well for
anisotropic particle in viscous shear flow with Navier slip boundary condition
at the particle surface, showing that the boundary slip can effectively
enhance the anisotropy of the particle, which fairly matches with the results
of [19]. In the above sections we have been assuming that $\Omega$ contains
only one particle. However, we can easily consider cases where $\Omega$
contains a large number particles. In the future, we will generalize our
method to the simulation of multiphase flow and particles interactions with
boundary slip.

## Acknowledgments

The authors would like to thank Prof. Roland Glowinski and Prof. Xiao-Ping
Wang for their useful discussions. This research is supported part by National
Key R & D Program of China (2018YFC0830300) and the National Natural Science
Foundation of China (No.11971020).

## References

## References

  * Feng et al. [1994] J. Feng, H. Hu, D. Joseph, Direct simulation of initial value problems for the motion of solid bodies in a Newtonian fluid, part 1. Sedimentation, J. Fluid Mech. 261 (1994) 95–134.
  * Glowinski [2003] R. Glowinski, Finite element methods for incompressible viscous flow in Handbook of Numerical Analysis, Vol. IX, North-Holland, Amsterdam, 2003.
  * Glowinski et al. [2008] R. Glowinski, J. He, J.L. Lions, Exact and approximate controllability for distributed parameter systems: A numerical approach, Cambridge University Press, Cambridge, UK, 2008.
  * Glowinski and He [2011] R. Glowinski, Q. He, A least-squares/fictitious domain method for linear elliptic problems with Robin boundary conditions, Commun. Comput. Phys. 9(3) (2011) 587–606.
  * Glowinski et al. [1994a] R. Glowinski, T. Pan, J.Periaux, A fictitious domain mehod for Dirichlet problems and applications, Comput. Methods Appl. Mech. Engrg. 111 (1994a) 283–303.
  * Glowinski et al. [1994b] R. Glowinski, T. Pan, J.Periaux, A fictitious domain method for external incompressible flow modeled by Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg. 112 (1994b) 133–148.
  * Glowinski et al. [1995a] R. Glowinski, T. Pan, J.Periaux, A Lagrange multiplier/fictitious domain method for the Dirichlet problem. Generalization to some flow problems, Japan J. Indust. Appl. Math. 12 (1995a) 87–108.
  * Glowinski et al. [1995b] R. Glowinski, T. Pan, A. Kearsley, J.Periaux, Numerical simulation and optimal shape for viscous flow by a fictitious domain method, Int. J. Numer. Methods Fluids 20 (1995b) 695–711.
  * He et al. [2011] Q. He, R. Glowinski, X.P. Wang, A least-squares/finite element method for the numerical solution of the Navier-Stokes-Cahn-Hilliard system modeling the motion of the contact line, J. Comput. Phys. 230 (2011) 4991–5009.
  * He et al. [2018] Q. He, R. Glowinski, X.P. Wang, The least-squares/fictitious domain method for incompressible viscous flow around obstacles with Navier slip boundary condition, J. Comput. Phys. 366 (2018) 281–297.
  * HH. [1996] H. HH., Direc simulation of flows of solid-liquid mixtures, Int. J. Multiphase Flow. 22 (1996) 335–352.
  * Jeffery [1922] G.B. Jeffery, The motion of ellipsoidal particles immersed in a viscous fluid, in Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, The Royal Society 102 (1922) 161–179.
  * MJ. [1992] H.H.J.D.C. MJ., Direc simulation of fluid particle motion, Theor. Comp. Fluid Dyn. 3 (1992) 285–306.
  * Pan et al. [1999] R.G.T. Pan, T. Hesla, D.D. Joseph, A distributed lagrange multiplier/fictitious method for flows around moving rigid bodies: application to particulate flows, Int. J. Multiphase Flow. 25 (1999) 755–794.
  * Pan et al. [2002] T.W. Pan, R. Glowinski, G. Galdi, Direct simulation of the motion of settling ellipsoid in Newtonian fluid, J. Comput. Appl. Math. 149 (2002) 71–82.
  * Peskin [1977] C. Peskin, Numerical ananlysis of blood flow in the heart, J. Comput. Phys. 25 (1977) 220–252.
  * Uhlmann [2005] M. Uhlmann, An immersed boundary method with direct forcing for the simulation of particulate flows, J. Comput. Phys. 209(2) (2005) 448–476.
  * Wang et al. [2020] Z. Wang, H. Q. He, J. Huang, The immersed boundary–lattice Boltzmann method for solving solid–fluid interaction problem with Navier–slip boundary condition (2020) submitted.
  * Zhang et al. [2015] J. Zhang, X. Xu, T. Qian, Anisotropic particle in viscous shear flow: Navier slip, reciprocal symmetry, and Jeffery orbit, Phys. Rev. E 91(3) (2015) 033016.