%\documentclass[lineno]{JFM-FLM_Au}
\documentclass[a4paper,12pt]{article}

% Fonts for LuaLaTeX
%\usepackage{fontspec}


% Math packages (amsmath first)
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsfonts}

% Symbols
\usepackage{wasysym}

% Tables
\usepackage{multirow}
\usepackage{tabularx}
\usepackage{adjustbox}

% Graphics and TikZ
\usepackage{import}
\usepackage{tikz, tikz-3dplot}
\usetikzlibrary{shapes,arrows.meta,chains,shapes.multipart,spy,backgrounds}
\usepackage{pgfplots}
\usepgfplotslibrary{fillbetween} 
\pgfplotsset{compat=1.18}
\usepackage{pgfplotstable}
\usepackage{xcolor}
\usepackage{overpic}
\usepackage{subcaption}
\usepackage{float}

% Lists
\usepackage{enumitem}

% Comments
\usepackage{comment}

% Bibliography
\usepackage{doi}
\usepackage{url}
\usepackage[numbers,sort&compress]{natbib}
\bibliographystyle{plainnat}

% Author/Affiliation
\usepackage{authblk}

\newcommand{\Pe}{\mathcal{P}}
\newcommand{\C}{\mathcal{C}}
\newcommand{\erfc}{\textrm{erfc}}
\newcommand{\ksi}{\xi}
\newcommand{\uz}{u^\star}
\newcommand{\ur}{u}
\newcommand{\GF}[1]{\mathcal{#1}}
\newcommand{\xio}{(\vec{\xi})}
\newcommand{\OKL}[1]{\left(#1\right)}
\newcommand{\ra}[1]{\renewcommand{\arraystretch}{#1}}

\newcommand{\red}[1]{{\color{red}#1}}
\newcommand{\tikzcircle}[2][red,fill=red]{\tikz[baseline=-0.5ex]\draw[#1,radius=#2] (0,0) circle ;}%
\let\oldRe\Re
\renewcommand{\Re}{\ensuremath{\text{Re}}}
\newcommand{\Ga}{\ensuremath{\text{Ga}}}
\def\gz #1{\mbox{\boldmath $\mit #1$}}


%%% TO AVOID CLASH BETWEEN SPY AND FILLBETWEEN LIBRARIES
\tikzset{spy scope/.append code={
   \patchcmd\tikz@atend@scope{\egroup\egroup}{%
      \pgf@remember@layerlist@globally%
      \egroup\egroup%
      \pgf@restore@layerlist@from@global%
   }{}{}
}}

\newtheorem{lemma}{Lemma}
\newtheorem{corollary}{Corollary}


\title{\vspace{-2em} % move up a bit if needed
Fluctuations, Clustering, and Interaction-Driven Dynamics in Sedimenting Particles at Low Galileo Numbers: A Neural Network Approach
\\[1ex] % space between title and the line below
\large\textit{Submitted to the Journal of Fluid Mechanics (under review)}}

\author[1]{Nejc Vovk\thanks{Corresponding author: nejc.vovk@um.si}}
\author[2]{Jana Wedel}
\author[2,3]{Paul Steinmann}
\author[1]{Jure Ravnik}

\affil[1]{Faculty of Mechanical Engineering, University of Maribor, Slovenia}
\affil[2]{Institute of Applied Mechanics, University of Erlangen-Nürnberg, Germany}
\affil[3]{Glasgow Computational Engineering Center, University of Glasgow, Scotland}

\date{}

\begin{document}
\maketitle

\begin{abstract}
In this study, we investigate the behaviour of sedimenting solid particles and the influence of microscopic particle dynamics on the collective motion of a sedimenting cloud. Departing from conventional direct numerical simulations (DNS), we introduce a novel machine learning framework, the Interaction-Decomposed Neural Network (IDNN), to model hydrodynamic particle interactions. The IDNN acts as a black-box module within a Lagrangian solver, predicting the particle drag force based on the relative positions of the nearest neighbours. This enables the recovery of force fluctuations, capturing effects previously accessible only through DNS. Our results show an increase in collective settling velocity in the dilute regime, consistent with earlier experimental and numerical studies, which we attribute to (i) fluctuations in the streamwise particle force around a value that is lower than the Stokes limit and (ii) the formation of particle clusters sedimenting at enhanced velocities. These fluctuations originate from persistent entrainment and ejection of particles in and out of the long, diffusive wakes generated by upstream particles at low Galileo numbers. Energy spectra of particle velocity fluctuations reveal a scale-dependent transfer of fluctuation energy, analogous to a turbulent-like cascade, with pronounced large-scale fluctuations at higher volume fractions. At very low volume fractions, fluctuation intensity and energy spectrum amplitudes diminish, though hydrodynamic interactions still remain appreciable.

%In this paper, we make use of a machine learning surrogate to try to explain the complex behaviour that occurs during settling of solid particles. This work introduces the Interaction-Decomposed Neural Network (IDNN), a new physics-based data-driven framework for estimating hydrodynamic particle interactions, whose essence is: (i) a two-block architecture, with a $1^{\text{st}}$ order block for pairwise interactions and a $2^{\text{nd}}$ block for higher-order effects; and (ii) input sorting by radial distance, implicitly encoding the relative influence of each neighbour. The developed physics-based data-driven model provides insights into the underlying mechanisms of particle settling that were previously accessible only through direct numerical simulations. We integrate the IDNN into the Lagrangian point-particle solver, to simulate the settling of solid particles and investigate the mechanisms behind the observed increase in settling velocity. The results suggest that particle clustering may contribute to this enhancement, even at low Galileo numbers.
\end{abstract}

%\begin{keywords}
%Particle sedimentation, Multibody approach, Hydrodynamic particle interaction, Neural Networks
%\end{keywords}
\noindent\textbf{Keywords:} Particle sedimentation, Multibody approach, Hydrodynamic particle interaction, Neural Networks

%%%%%%%%%%%%%%%%%%%%
\import{}{introduction.tex}
\import{}{objective_outline.tex}
\import{}{deterministic_model.tex}
\import{}{particleLadenFlow.tex}
\import{}{results.tex}
\import{}{sedimentation.tex}
\import{}{conclusions.tex}
%%%%%%%%%%%%%%%%%%%%
\clearpage

\section*{Appendix}
\appendix
%\begin{comment}
%\begin{appen}
\section{Particle force computation with BEM}\label{appA}
The governing equation for the steady, incompressible flow of a Newtonian fluid is solved, as described in Eq. (\ref{eq:StokesEquation}). The Stokes flow Green's functions satisfy the continuity equation $\mathbf{\nabla}\cdot\mathbf{u}_\text{f} = 0$ and are the solutions of the singularly forced Stokes equation.
The 3D free-space Green's functions are
\begin{equation}\label{e:StkGF}
	\GF{G}^\star_{ij}=\frac{\delta_{ij}}{r}+\frac{\hat r_i\hat r_j}{r^3}, \qquad
	\GF{T}^\star_{ijk}=-6\frac{\hat r_i\hat r_j\hat r_k}{r^5}.
\end{equation}
The boundary integral representation for the Stokes problem is \citep{pozrikidisIntroductionTheoreticalComputational2011}:
\begin{equation}\label{eq5656}
	c(\gz \ksi) u_j(\gz \ksi) = \int_{\Gamma}^{PV}u_i\GF{T}^\star_{ijk}n_k d\Gamma
	-\frac{1}{\mu}\int_\Gamma \GF{G}^\star_{ji} q_i \text{d}\Gamma,
\end{equation}
where $c(\gz \ksi)=2\alpha$ is twice the solid angle as seen from the point $\gz \ksi$, i.e. in the interior of the domain $c=8\pi$, at a smooth boundary $c=4\pi$. The boundary tractions are denoted by $\mathbf{q} = \gz \sigma\cdot \mathbf{n}$. The normal vector $\mathbf{n}$ points into the domain. The terms on the right represent the double and single layer potentials of the three-dimensional Stokes flow.
To derive a discrete version of (\ref{eq5656}) we consider the boundary $\Gamma = \sum_l\Gamma_l$ to be decomposed into boundary elements $\Gamma_l$:
\begin{equation}
	c(\gz \ksi) u_j(\gz \ksi) = 
	\sum_l\int_{\Gamma_l}^{PV}u_i\GF{T}^\star_{ijk}n_k^{(l)} \text{d}\Gamma
	-\frac{1}{\mu}\sum_l\int_{\Gamma_l} \GF{G}^\star_{ji}q_i \text{d}\Gamma,
\end{equation}
where $n_k^{(l)}$ is the $k$ component of the normal vector pointing from boundary element $l$ into the domain.

Let $\Phi$ be the interpolation functions used to interpolate the function values within boundary elements, i.e. $u_i=\sum_m\Phi_m u_i^{(l,m)}$, where $u_i^{(l,m)}$ is the $m^{th}$ nodal value of function within $l^{th}$ boundary element. Constant interpolation is considered for flux. This yields:
\begin{eqnarray}
	c (\gz \ksi) u_j (\gz \ksi) = 
	\sum_l\sum_mu_i^{(l,m)}\int_{\Gamma_l}^{PV}\Phi_m\GF{T}^\star_{ijk}n_k^{(l)}\text{d}\Gamma %\nonumber \\
	-\frac{1}{\mu}\sum_lq_i^{(l)}\int_{\Gamma_l} \GF{G}^\star_{ji} \text{d}\Gamma.
\end{eqnarray}
The following integrals must be calculated for each boundary element $l$:
\begin{eqnarray}\label{eq_int}
	T_{ij}^{(l,m)}(\mathbf{\ksi}) = \int_{\Gamma_l}^{PV} \Phi_m \GF{T}^\star_{ijk}n_k^{(l)}\text{d}\Gamma,
	\nonumber \hspace{0.5cm}
	G_{ij}^{(l)}(\gz \ksi) = \int_{\Gamma_l} \GF{G}^\star_{ij} \text{d}\Gamma.
\end{eqnarray}
Considering boundary conditions we can place the source point into nodes, where unknown values are located and produce a system of linear equations for the velocity and traction. The Andromeda code is able to efficiently simulate Stokes flow based on boundary only discretization. As such it is ideally suitable for performing numerous simulations needed to develop ML based models, as is the subject of present research.
 
Computationally the most expensive part of the simulation is finding the solution of the system of linear equations, created by the BEM based discretization procedure. To facilitate the possibility of parallel computing, we use the  {\it mpich} library  to set up the system of linear equations in parallel and the {\it LIS} library \citep{nishidaExperienceDevelopingOpen2010} to find the solution also in parallel.

\section{Mesh validation study}\label{appC}
For this analysis, we focus on the discretization of a single particle in a plug flow and compare the simulated drag force with the analytical solution of the Stokes drag, Eq. (\ref{eq:stokesDrag}). The computational domain is identical to that shown in Fig. \ref{fig:problemDefinition}. A Dirichlet boundary condition is applied to the velocity field on the outer sphere to simulate plug flow, and on the surface of the particle to enforce a no-slip condition. A Neumann boundary condition is imposed on the particle surface for the pressure field. The results plotted in Fig. \ref{fig:plugFlowMeshValidation} show good convergence and the chosen mesh density satisfies both the conditions of good accuracy and computational affordability. For subsequent simulations, where more than one particle is considered in the flow, we keep the mesh design for all particles the same as the particle mesh in the validation study, marked in red. This domain mesh along with the discretized particle, is shown in Fig.~\ref{fig:finalMesh}. We further quantitatively assess the discretization uncertainty by using the method proposed by \citet{celikProcedureEstimationReporting2008}. The BEM numerical method expresses a strong monotone convergence of order $p=2.52$. The numerical uncertainty, in terms of the grid convergence index (GCI), accounts to $8.18\%$. Detailed results are presented in Tab. \ref{tab:GCIresults}. Since the mesh for each of the simulations changes due to changing particle positions, we automized the meshing procedure via Python scripts calling the {\it gmsh} \citep{geuzaineGmsh3DFinite2009} mesher. 
\begin{figure}[h]
    \centering
    \begin{tikzpicture}[scale=0.75]
        \begin{axis}
            [
                ylabel={$\text{Re}_\text{p} c_\text{D}$ (-)},
                xlabel={Number of nodes},
                grid=major,
                width = 0.35\textwidth,
                xmin =1,
                xmax=1700,
                ymax=25,
                legend style={at={(1.05,1)}, anchor=north west},        xticklabel style={
                    /pgf/number format/.cd,
                    1000 sep={}  % thin space for thousands separator
                },
                scaled x ticks=false 
            ]

            \addplot[
                only marks, 
                mark=o, 
                color=black
            ] 
            table [col sep=semicolon] {plots/streamwiseForce.csv};
            \addlegendentry{BEM simulation}

            \node[label={\footnotesize $88\%$},circle,fill,inner sep=0.5pt] at (axis cs:212.04141747746138, 10.043483250659897) {};
            \node[label=south:{\footnotesize $92\%$},circle,fill,inner sep=0.5pt] at (axis cs:346.16365514339276, 17.522500350671372) {};
            \node[label=south east:{\footnotesize $78\%$},circle,fill,inner sep=0.5pt] at (axis cs:547.9144616875583, 19.818875046224864) {};
            \node[label=west:{\footnotesize $77\%$},circle,fill,inner sep=0.5pt] at (axis cs:827.1126356460643, 22.189770597161473) {};
            \node[label=south:{\footnotesize $70\%$},circle,fill,inner sep=0.5pt] at (axis cs:1088.3691868249578, 23.245304191479324) {};
            \node[label=south:{\footnotesize $72\%$},circle,fill,inner sep=0.5pt] at (axis cs:1506.7201387383484, 24.380867369709648) {};

            \addplot[
                color=gray,
                style=dashed
            ] table[col sep=space, header=true] {
                X    Y  
                0    24
                2000 24    
            };
            \addlegendentry{Stokes drag}

            \addplot[
                only marks, 
                mark=*,
                mark size=3pt,
                color=red
            ] coordinates {(827.1126356460643, 22.189770597161473)};

        \end{axis}
    \end{tikzpicture}
        \caption{Force exerted by the fluid on a single particle during plug flow versus the number of mesh nodes used. The symbol labels refer to the share of nodes used to discretize the particle, while the rest was used to discretize the outer spherical domain. The mesh chosen for further simulations is shown in red.}
        \label{fig:plugFlowMeshValidation}
\end{figure}
\begin{table}
    \begin{center}
  \def~{\hphantom{0}}
    \begin{tabular}{llr}
        \textbf{Parameter} & \textbf{Symbol} & \textbf{Value} \\[3pt]
        Order of convergence & $p$ & $2.52$ \\
        Coarse mesh extrapolated result & $(\text{Re}_\text{p} c_\text{D})_{\text{ext}, 32}$ & $24.86$ \\
        Fine mesh extrapolated result & $(\text{Re}_\text{p} c_\text{D})_{\text{ext}, 21}$ & $26.32$ \\
        Coarse mesh numerical uncertainty & $\text{GCI}_{\text{coarse}, 32}$ & $2.95\%$ \\
        Fine mesh numerical uncertainty & $\text{GCI}_{\text{fine}, 21}$ & $8.18\%$
    \end{tabular}
    \caption{Results of GCI analysis for plug flow over a single particle.}
    \label{tab:GCIresults}
    \end{center}
\end{table}  
\begin{figure}
    \centering
    \includegraphics{figures/final-mesh.pdf}
    \caption{The mesh, recognized as a good compromise between the accuracy and the computational cost, that was used for running numerous simulations during the training database generation. The colour on the particle surface demonstrates the pressure distribution on the particle surface, as a result of the BEM simulation.}\label{fig:finalMesh}
\end{figure}

\section{Coordinate system transformations}\label{appB}
We observe a cloud of $N$ particles in the fluid flow, where each reference particle, denoted as $i$, is surrounded by a cluster of $M$ closest neighbours, denoted as $j$. The short inter--particle distance causes interactions of the surrounding flow fields, resulting in a disturbance of the reference particle drag force. We consider the cloud of particles in two coordinate system definitions. The global coordinate system (GCS) denotes the global coordinates of the reference particle, 
\begin{equation}
    \underline{r}_i = \underline{e}_1 x_i + \underline{e}_2 y_i + \underline{e}_3 z_i,     
\end{equation}
and its neighbours 
\begin{equation}
    \underline{r}_{i,j} = \underline{e}_1 x_{i,j} + \underline{e}_2 y_{i,j} + \underline{e}_3 z_{i,j}, 
\end{equation}
where $\underline{e}_1 \dots \underline{e}_3$ form the orthonormal base of the GCS and are defined as $\underline{e}_1 = [1, 0, 0]$, $\underline{e}_2 = [0, 1, 0]$ and $\underline{e}_3 = [0, 0, 1]$. The second considered coordinate system is the local coordinate system (LCS) of the reference particle, with the corresponding coordinates denoted as 
\begin{equation}
    \underline{r}^{\hspace{2pt} '}_i = [0, 0, 0],
\end{equation}
\begin{equation}
    \underline{r}^{\hspace{2pt} '}_{i,j} = \underline{e}^{\hspace{2pt} '}_{1, i} x^{\hspace{2pt} '}_{i,j} + \underline{e}^{\hspace{2pt} '}_{2, i} y^{\hspace{2pt} '}_{i,j} + \underline{e}^{\hspace{2pt} '}_{3, i} z^{\hspace{2pt} '}_{i,j},
\end{equation}
where $\underline{e}^{\hspace{2pt} '}_{1, i} \dots \underline{e}^{\hspace{2pt} '}_{3, i}$ form the orthonormal base for the LCS for each reference particle. The subscripts $i,j$ in the above definitions denote that the coordinate corresponds to the $j$--th neighbour of the $i$--th reference particle. The reason behind considering two coordinate systems is that the whole training dataset is defined in the LCS, where the base vector $\underline{e}^{\hspace{2pt} '}_{2, i}$ is aligned with the relative velocity vector at the position of the reference particle, $\underline{e}^{\hspace{2pt} '}_{2, i} || \underline{u}_{rel, i}$, as shown in Fig. \ref{fig:GCS-LCS}.
\begin{figure}[h]
    \centering
    \begin{tikzpicture}
        \draw[-Stealth] (0,0) -- (1,0) node[below left] {$\underline{e}_1$};
        \draw[-Stealth] (0,0) -- (0,1) node[below left] {$\underline{e}_2$};
        \draw[-Stealth] (0,0) -- (-135:0.707) node[below left] {$\underline{e}_3$};

        \draw[-Stealth, dashed] (3,1) -- (1,2) node[right, xshift=10pt, yshift=-1pt] {$\underline{r}_{i,j}^{\hspace{2pt} '}$};
        \draw[-Stealth, dashed] (0,0) -- (1,2) node[left, xshift=-8pt] {$\underline{r}_{i,j}$};
        \draw[ball color=cyan!50!blue, opacity=0.5, draw opacity=1] (1,2) circle[radius={0.3cm}];
        \draw[-Stealth, dashed] (0,0) -- (3,1) node[left, xshift=-15pt] {$\underline{r}_{i}$};

        \draw[ball color=cyan!50!blue, opacity=0.5, draw opacity=1] (3,1) circle[radius={0.3cm}];
        \draw[-Stealth] (3,1) -- ++ (-10:2) node[right] {$\underline{u}_{rel, i}$};
        \draw[-Stealth] (3,1) -- ++ (-10:1) node[above] {$\underline{e}_{2,i}^{\hspace{2pt} '}$};
        \draw[-Stealth] (3,1) -- ++ (-100:1) node[right] {$\underline{e}_{1,i}^{\hspace{2pt} '}$};
        \draw[-Stealth] (3,1) -- ++ (100:0.707) node[above] {$\underline{e}_{3,i}^{\hspace{2pt} '}$};        
    \end{tikzpicture}
    \caption{Visualization of the GCS and the LCS.}
    \label{fig:GCS-LCS}
\end{figure}

It can be seen, that in order to obtain the reference particle force in the global coordinate system, a series of transformations has to be applied to the global neighbour particle coordinates. The transformation of the neighbour particle position vector from GCS to LCS can be generally written as
\begin{equation}\label{eq:GCS-LCS}
    \underline{r}^{\hspace{2pt} '}_{i,j} = \underbar{R}_i \left[ \underline{r}_{i,j} - \underline{r}_i \right],
\end{equation}
where $\underbar{R}_i$ is the rotation matrix for the $i$--th reference particle. The rotation matrix has to be constructed, so that the collinearity between the relative velocity vector of the reference particle and the base vector $\underline{e}^{\hspace{2pt} '}_{2, i}$ is satisfied. The rotation matrix can be constructed for two linearly independent vectors, using the Rodrigues' rotation formula \citep{wellerTensorialApproachComputational1998}, which in our case reads as
\begin{equation}
    \underbar{R}_i = 
        c \underbar{I} + \left[ \frac{\underline{u}_{rel, i}}{| \underline{u}_{rel, i} |} \otimes \underline{e}^{\hspace{2pt} '}_{2, i} - 
        \underline{e}^{\hspace{2pt} '}_{2, i} \otimes \frac{\underline{u}_{rel, i}}{| \underline{u}_{rel, i} |} \right] + 
        \left[ 1 - c \right] \frac{ \underline{a} \otimes \underline{a} }{|\underline{a}|^2 },
\end{equation}
where 
\begin{equation}
    c = \underline{e}^{\hspace{2pt} '}_{2, i} \cdot \frac{\underline{u}_{rel, i}}{| \underline{u}_{rel, i} |}
\end{equation}
and
\begin{equation}
    \underline{a} = \underline{e}^{\hspace{2pt} '}_{2, i} \times \frac{\underline{u}_{rel, i}}{| \underline{u}_{rel, i} |}.
\end{equation}
In above equations, the operators $\otimes$, $\cdot$ and $\times$ represent the dyadic product, dot product and the cross product respectively. The above rotation matrix definition holds if $\underline{e}^{\hspace{2pt} '}_{2, i}$ and $\frac{\underline{u}_{rel, i}}{| \underline{u}_{rel, i} |}$ are linearly independent. If the vectors are collinear and contradirectional ($c < 0$), the rotation matrix is constructed as
\begin{equation}
    \underbar{R}_i = - \underbar{I} + 2 \frac{\underline{b} \otimes \underline{b}}{|\underline{b}|},
\end{equation}
where $\underline{b}$ is a vector, perpendicular to $\underline{e}^{\hspace{2pt} '}_{2, i}$. In the case where $\underline{e}^{\hspace{2pt} '}_{2, i}$ and $\frac{\underline{u}_{rel, i}}{| \underline{u}_{rel, i} |}$ are collinear and codirectional ($c > 0$), the rotation matrix is equal to the identity,
\begin{equation}
    \underbar{R}_i = \underbar{I}.
\end{equation}
To be able to use the obtained force prediction in the Lagrangian solver, the obtained prediction must be transformed with the rotation matrix back to the GCS as
\begin{equation}
    \underline{F}_{\text{IDNN}, i} = \underline{R}_i^\top \underline{F}'_{\text{IDNN}, i}.
\end{equation}
%\end{appen}
%\end{comment}
\clearpage

\textbf{Acknowledgements.} The authors would like to thank the Slovenian Research and Innovation Agency (research core funding No. P2-0196 and project J7-60118) and the Deutsche Forschungsgemeinschaft (project STE 544/75-1).

\textbf{Declaration of Interests.} The authors report no conflict of interest.

%\bibliographystyle{plain} 
\bibliography{/home/nejcv/Documents/PROJEKTI/bibliography/library}

\end{document}
