A finite element / neural network framework for modeling suspensions of non-spherical particles. Concepts and medical applications

An accurate prediction of the translational and rotational motion of particles suspended in a fluid is only possible if a complete set of correlations for the force coefficients of fluid-particle interaction is known. The present study is thus devoted to the derivation and validation of a new framework to determine the drag, lift, rotational and pitching torque coefficients for different non-spherical particles in a fluid flow. The motivation for the study arises from medical applications, where particles may have an arbitrary and complex shape. Here, it is usually not possible to derive accurate analytical models for predicting the different hydrodynamic forces. However, considering for example the various components of blood, their shape takes an important role in controlling several body functions such as control of blood viscosity or coagulation. Therefore, the presented model is designed to be applicable to a broad range of shapes. Another important feature of the suspensions occurring in medical and biological applications is the high number of particles. The modelling approach we propose can be efficiently used for simulations of solid-liquid suspensions with numerous particles. Based on resolved numerical simulations of prototypical particles we generate data to train a neural network which allows us to quickly estimate the hydrodynamic forces experienced by a specific particle immersed in a fluid.


Introduction
The prediction of the motion of non-spherical particles suspended in a fluid is crucial for the understanding of natural processes and industrial applications. In such processes, particles can have different shapes and sizes, may be deformed and can interact with each other. So far, in the majority of scientific studies, particulate flow modelling is investigated with the hypothesis of perfectly spherical particles, thereby eliminating orientation and shape effects. This assumption is very convenient due to its simplicity, the fact that the behaviour of spheres is well known and the availability of a number of models to describe the interaction with fluid flow. The study of suspensions of multiple, irregular-shaped, interacting and deformable particles has received less attention and still presents a challenge.
Particles come in all sort of shapes and sizes, in fact, due to the arbitrary nature of naturally occurring particles there are an indefinite number of possible shapes. On the other hand, there is a common understanding that particle shape has a strong influence on the dynamics of NSPS (non-spherical particulate systems). These two factors combined makes modelling of NSPS in general way impossible, since for describing the motion of non-spherical particles, detailed information on the fluid dynamic forces acting on such particles are necessary, but generally not available. Therefore, particular models place emphasis on different shapes and types of flow. In our work we focus on medical applications, namely on modelling of platelets dynamics under blood flow. Platelets play a main role in the process of blood coagulation and therefore are of great interest in the modelling of blood flow. The majority of models characterize platelet motion quantitatively and use approaches such as immersed boundary method [13], cellular Potts model [39,40,41] and dissipative particle dynamics [12,35], treating platelets as points and thus neglecting entirely their shape. Some effort has been done to model platelets as rigid twoor three-dimensional spheres or spheroids [23,34,45] but this simplification of shape has been shown to affect processes in which platelets are involved (e.g. spherical platelets marginate faster than ellipsoidal and disc forms) [36].
On the other hand, several studies have been performed to model particles of irregular shapes, but motivation behind them usually arises from engineering applications (dispersion of pollutant, pulverized coal combustion, pneumatic transport) [46], where particles are much bigger and usually constitute a significant part of the volume of the suspension [33]. But even under given very specific circumstances there is no set of correlations of the forces acting on irregular-shaped particles suspended in a fluid (forces arising from fluid-particle interaction). Furthermore, in these kinds of models interactions between particles become dominant in terms of determining particle dynamics, whereas platelets are very dilute in blood and their contact is rather rare.
The motivation behind this work is fourfold and arises from the specificity of the modelled phenomenon. Firstly, platelets, because of their small size (compared to blood vessel), are often modelled quantitatively, ignoring the importance of shape effects. Secondly, platelets constitute only a very small volume of blood what makes most of the engineering applications-driven models, where particles are very dense and often interact with each other, inappropriate. Thirdly, highly irregular shape of platelets requires new methods for estimating force coefficients in fluid-particle interaction. Fourthly, platelets are numerous and the evolution of their movement needs to be evaluated quickly and efficiently.
A common approach to the problem of modelling of NSPS consists in developing analytical models for fluid dynamic forces acting on particles, cf. [16,20,22,44,46]. In this contribution we go a different way and refrain from giving analytical expressions. Instead, by prototypical simulations we train a neural network model that takes a couple of parameters describing the particle shape, size and the flow configuration as input and which gives hydrodynamical coefficients like drag, lift and torque as output. We demonstrate that such a model can be efficiently trained in an offline phase for a range of particles. Lateron, the coupling of the flow model with the particle system only requires the evaluation of the network for getting updates on these coefficients. This two-step approach with an offline phase for training a network based on the particle classes under consideration allows for a direct extension to further applications.
A similar approach is considered by in [42]. Here, the authors design and train a radial basis function network to predict the drag coefficient of nonspherical particles in fluidized beds. Training is based on experimental data and the network input is the particle's sphericity and the Reynolds number, covering the Stokes and the intermediate regime. Our approach, based on training data generated by detailed simulations of prototypical particles for predicting drag, lift and torque coefficients could by augmented by including further experimental data.
In the following section we will describe prototypical medical applications where such a heterogenous modeling approach can be applied. Then, in Section 3 we detail the general framework for coupling the Navier-Stokes equation with a discrete particle model. Section 4 introduces the neural network approach for estimating the hydrodynamical coefficients and we describe the procedure for offline training of the network. Then, different numerical test cases are described in Section 5.1 that show the potential of such a heterogeneous modeling approach. Section 5.2 is devoted to a numerical study comprising many particles and shows the efficiency of the presented approach. We summarize with a short conclusion in Section 6.

Modeling of suspensions with non-spherical particles and medical applications
Platelets are a vital component of the blood clotting mechanism. They are small non-nucleated cell fragments. They have a diameter of approximately 2 − 4 µm, thickness of 0.5 µm, volume of about 7 µm 3 and a number density of 1.5 − 4 · 10 5 µl −1 [11] which leads to a volume fraction of only about 10 −3 : 1.
Still they are a vital component of the blood clotting mechanism. In the rest state platelets shape is discoid, but they have the ability of deforming as a response to various stimuli (chemical and mechanical). They may become star shaped (rolling over blood vessels wall to inspect its integrity). During the clotting process they undergo deep morphological changes, from becoming spherical and emitting protuberances (philopodia or pseudopods) which favour mutual aggregation, as well as adhesion to other elements constituting the clot to fully flat, spread stage, to enable wound closure. Thrombocytes constitute approximately less than 1% of the blood volume, therefore individual platelets have a negligible effect on blood rheology [19].
Due to a significant effect of the particles shape on their motion and their practical importance in industrial applications, the non-sphericity started attracting attention in the modelling and simulations of particles transport in fluid flows [17,22,44]. Unfortunately, it is not possible to consider each shape in the implementation of numerical methods because of non-existence of a single approach describing accurately the sizes and shapes of non-spherical particles. Spheres can be described by a single characteristic value, i.e. the diameter, whereas non-spherical particles require more parameters. Even very regular shapes need at least two parameters. Moreover, the particles may have varying orientation with respect to the flow, what makes the description of their behaviour even more difficult. Even though several methods for shape parametrization and measurement have been suggested, none has won greater acceptance. One of most commonly used shape factor is a sphericity which was firstly introduced in [37] and defined as the ratio between the surface area of a sphere with equal volume as the particle and the surface area of the considered particle. Then the drag coefficient for non-spherical particle is estimated by using correlations for spherical particles and modified to take into account the sphericity factor [43]. Models using sphericity as a shape factor give promising results when restricted to non-spherical particles with aspect ratios β less than 1.7 [7], where β = L/D with L and D being length and diameter of the considered particle. For particles having extreme shapes and those having little resemblance to a sphere, the sphericity concept fails to produce satisfying quantitative results [6]. In general, the lower the sphericity, the poorer is the prediction. Also, the same value of the sphericity may be obtained for two very varying shapes whose behaviour in the flow is different. Moreover, the sphericity does not take the orientation into account. In order to introduce orientation dependency in drag correlations, some researchers use two additional factors: the crosswise sphericity and lengthwise sphericity [16]. Most of these correlations employ also dependency on particle Reynolds number defined as where ρ and µ are fluid density and viscosity,ū = u f − u p is the velocity of the particle relative to the fluid velocity and d eq is the equivalent particle diameter, i.e. the diameter of a sphere with the same volume as the considered particle in order to include the importance of fluid properties.
The shape factor concept may be described as an attempt to define a single correlation for drag for all shapes and orientations. Another approach appeared as an alternative consisting in obtaining drag coefficient expressions for a fixed shape and any orientations: the drag coefficient is determined at two extreme angles of incidence (0 • and 90 • ) from existing correlations which are then linked by some functions resulting in the whole range of angles of incidence for non-spherical particles [30]. However, besides drag force, nonspherical particles are associated with orientation and shape induced lift along with pitching and rotational torques. Hölzer and Sommerfeld [17] investigated a few different shapes of non-spherical particles at different flow incident angles using the Lattice Boltzmann method to simulate the flow around the particle. Wachem et al. [44] proposed new force correlations (for drag, lift, pitching and rotational torque) for particular shapes of non-spherical particles (two ellipsoids with different aspect ratio, disc and fibre) from data given by a direct numerical simulation (DNS) carried out with an immersed boundary method. Those correlations employ particle Reynolds number, angle of incidence and some shape-related coefficients. Ouchene et al. [25] determined force coefficients depending on particle Reynolds number, aspect ratio and angle of incidence by fitting the results extracted from DNS computations of the flow around prolate ellipsoidal particles. Discrete element methods (DEM) coupled with computational fluid dynamics (CFD) has been recognized as a promising method to meet the challenges of modelling of NSPS [47,48]. DEM is a numerical approach for modelling a large number of particles interacting with each other. The simplest computational sequence for the DEM typically proceeds by solving the equations of motion, while updating contact force histories as a consequence of contacts between different discrete elements and/or resulting from contacts with model boundaries. It is designed to deal with very dense suspensions, where contacts between particles are very common and play a key role in determining the motion of particles, see [46] for an extensive overview of DEM. Obviously, for a particle with a certain specific shape, the general expressions derived from the first factor shape approach tend to be less accurate than the specialized one for that shape, but the efficiency of interpolations/extrapolations to the various shapes to provide the general expression is an attractive perspective on engineering applications. On the other hand, particles occurring in biological processes are usually very numerous. Therefore, an effective method not only has to be accurate but also efficient in terms of computational time.
To overcome the aforementioned limitations in modelling of NSPS we employ the recently trending approach and use machine learning to design a method that enables us to model the behaviour of suspensions of particles of an arbitrary shape while maintaining at the same time the accuracy of shape-specific models. We also place an emphasis on the computational efficiency as usually there are plenty of particles involved in medical processes and engineering problems.

Model description
This section describes the general numerical framework for suspensions of particles in a Navier-Stokes fluid. The discretization of the Navier-Stokes equations is realized in the finite element toolbox Gascoigne 3D and outlined in Section 3.1. Then, in Section 3.2 we describe a very simple model for the motion of the particles.

Fluid dynamics
Consider a finite time interval I = [0, T ] and a bounded domain Ω ∈ R d for d ∈ {2, 3}. We assume incompressibility of fluid, which is modelled by the Navier-Stokes equations that take the form where v denotes the fluid velocity, σ is the Cauchy stress tensor p the pressure, ρ the fluid mass density and ν the kinematic viscosity. Fluid density and viscosity are assumed to be nonnegative and constant. The fluid boundary is split into an inflow boundary Γ in , an outflow boundary Γ out and rigid no-slip wall boundaries Γ wall . On the inflow and walls we impose Dirichlet boundary conditions while on the outflow we apply the donothing condition (see e.g. [15]) where v in is prescribed inflow-profile and n is the outward unit normal vector.
Discretization For temporal discretization of the Navier-Stokes equations we introduce a uniform partitioning of the interval I = [0, T ] into discrete steps 0 = t 0 < t 1 < · · · < t N = T, k := t n − t n−1 .
By v n := v(t n ) and p n := p(t n ) we denote the approximations at time t n . We use a shifted version of the Crank-Nicolson time discretization scheme which is second order accurate and which has preferable smoothing properties as compared to the standard version, see [14, Remark 1], i.e.
where, typically, θ = 1+k 2 . For spatial discretization we denote by Ω h a quadrilateral (or hexahedral) finite element mesh of the domain Ω that satisfies the usual regularity assumptions required for robust interpolation estimates, see [29,Section 4.2]. Adaptive meshes are realized by introducing at most one hanging node per edge. Discretization is based on second order finite elements for pressure and velocity. To cope with the lacking inf-sup stability of this equal order finite element pair we stabilize with the local projection method [2]. Local projection terms are also added to stabilize dominating transport [3]. Finally, velocity v n ∈ [V h ] 2 and pressure p n ∈ V h (where we denote by V h the space of bi-quadratic finite elements on the quadrilateral mesh) are given as solution to where the stabilization parameters are element-wise chosen as [5] we denote the interpolation into the space of bi-linear elements on the same mesh Ω h .
Solution of the discretized problem Discretization by means of (4) gives rise to a large system of nonlinear algebraic equation which we approximate by a Newton scheme based on the analytical Jacobian of (4). The resulting linear systems are solved by a GMRES iteration (Generalized minimal residual method [31]), preconditioned by a geometric multigrid solver [1]. As smoother we employ a Vanka-type iteration based on the inversion of the submatrices belonging to each finite element cell. These local 27 × 27 (108 × 108 in 3d) matrices are inverted exactly. Essential parts of the complete solution framework are parallelized using OpenMP, see [10].

Particle dynamics
The particles suspended in the fluid are described as rigid bodies and their dynamics is driven by the hydrodynamical forces of the flow. Each particle P with center of mass x P , velocity V P and angular velocity Ω P is governed by Newton's law of motion where m P is the particle's mass, and J P its moment of inertia given by with the (uniform) particle density ρ P . n P is the unit normal vector on the particle boundary facing into the fluid.
A resolved simulation is out of bounds due to the large number of platelets and in particular due to the discrepancy in particle diameter (about 10 −6 m) versus vessel diameter (about 10 −3 m). Instead, we consider all platelets to be point-shaped and determine traction forces F(v, p; P ) and torque T(v, p; P ) based on previously trained neural networks. These coefficients will depend on the shape and the size of the particles but also on their relative orientation and motion in the velocity field of the blood. Since the relative velocities (blood vs. particles) are very small the interaction lies within the Stokes regime with a linear scaling in terms of the velocity. The deep neural network will predict coefficients for determining coefficients of drag C d , lift C l , pitching torque C p and rotational torque C r and the resulting forces exerted on each particle P are given by where P = (L x , L y , L z , α top , α bot ) describes the particle shape and where ψ P is the relative angle of attack which depends on the particle orientation but also on the relative velocity vector between blood velocity and particle trajectory, see Figure 2. The coefficent functions C d , C l , C p and C r will be trained based on detailed numerical simulations using random particles in random configurations. By (v −V p ) ⊥ we denote the flow vector in lift-direction, orthogonal to the main flow direction. In 3d configurations, two such lift coefficients must be trained. Here we will however only consider 2d simplifications with one drag and one lift direction.

An artificial neural network model for predicting hydrodynamical parameters
In this section we describe the neural network model for coupling the Navier-Stokes equation with a suspension of non-spherical particles. The different hydrodynamical coefficients will be taken from a neural network, which is trained in an offline phase. Training data is achieved by resolved Navier-Stokes simulations using prototypical particles with random parameters. The setting investigated in this work carries several special characteristics that differ from industrial applications.
-The particle density is very small -about 1.04 − 1.08 ·10 3 gl −1 and the particle and fluid densities are similar (the average density of whole blood for a human is about 1.06 · 10 3 gl −1 ). Blood contains about 200 000 -400 000 platelets per mm 3 summing up to less than 1% of the overall blood volume [38]. Hence we neglect all effects of the particles onto the fluid. This simplification is possible since we only model the platelets as rigid particles. Effects of the red blood cells, much larger and appearing in greater quantity, can be integrated by means of a non-Newtonian rheology. -The particle Reynolds numbers are very small (with order of magnitude about 10 −4 or less) such that we are locally in the Stokes regime. This is mainly due to the smallness of the platelets (diameter approximately 3 µm) and the small flow velocities at (bulk) Reynolds numbers ranging from 50 to 1 000 depending on the specific vessel under investigation. We focus on coronary vessels with a diameter around 2 mm and with Reynolds number about 200. -The platelets have a strongly non-spherical, disc-like shape. Their shape and size underlies a natural variation. Furthermore, under activation, the particles will take a spherical shape.
Instead of deriving analytical models for the transmission of forces from the fluid to the particles, we develop a neural network for the identification of drag, lift and torque coefficients based on several parameters describing the shape and the size of the platelets and the individual flow configuration.

Parametrization of the platelets
We model the platelets as variations of an ellipsoid with major axes L x ×L y ×L z with L x ≈ L z ≈ 3 µm and L y ≈ 0.5 µm. In y-direction upper α top and lower α bot semi-ellipsoids are modified to give them a more or less concave or convex shape. Alltogether, each particle is described by a set of 5 parameters P = (L x , L y , L z , α top , α bot ). The surface of the platelets is given as zero contour of the levelset function where we define We assume that all parameters L x , L y , L z , α top , α bot are normally distributed with means indicated above and with standard derivation 0.3 for the lengths L x , L y , L z and 0.4 for the shape parameters α top , α bot . We drop particles that exceed the bounds In Fig. 1 we show some typical shapes of the platelets. Next, we indicate mass, center of mass and moment of inertia for a parametrized particle P = (L x , L y , L z , α top , α bot ). Unless otherwise specified, all quantities are given in µm and g. The mass of a particle P is approximated by This approximation is based on a weighted one-point Gaussian quadrature rule. It is accurate with an error of at most 2% for all α ∈ [0.2, 2]. The center of mass for a particle P is given by \ Fig. 2 Left: typical configuration of one platelet (in red) within velocity field of the blood (blue arrows). The angle ϕ is the orientation of the platelet relative to its standard orientation. Right: the black arrow δv = v − V is the velocity acting on the platelet. By ψ := ∠(δv, ex) − ϕ we denote the effective angle of attack.
The moment of inertia in the x/y plane (the only axis of rotation that we will consider in the 2d simplification within this work) is given by which is accurate up to an error of at most 1% for all α ∈ [0.2, 2]. These coefficients are computed once for each particle and stored as additional parameters.
2d Simplification To start with, we apply a two dimensional simplification of the problem by assuming that the blood vessel is a layer of infinite depth (in z-direction) and that it holds v 3 = 0 for the blood velocity and V 3 = 0 for all particles. Further, given the symmetry of the particles w.r.t. rotation in the x/y-plane, no traction forces in z-direction will appear. Besides that, rotation is restricted to rotation around the x/y-plane. Hence, Ω = (0, 0, ω) is described by a scalar component. A complete particle is then described by where (L x , L y , L z , α top , α bot ) are the 5 shape parameters and, X is the (2d) position, ϕ the orientation w.r.t. the z-axis, V the (2d) velocity and Ω the angular velocity w.r.t. the z-axis rotation.
To describe the forces acting on the particle suspended in the Navier-Stokes fluid we denote by δv := v − V the effective velocity vector, i.e. the relative velocity that is acting on the platelets. By ψ := ∠(e x , v − V) − ϕ we denote the effective angle of attack which is the angle between relative velocity δv and the current orientation of the platelet, see Fig. 2 (right) and (6). It is computed as We denote by v(x i ) the Navier-Stokes velocity at node x i , by v(t i ) the unit-tangential vector at this node in counter-clockwise orientation. With v i , t i we denote the velocity contribution in tangential direction. Positive values are indicated in orange (at node x 2 ), negative contributions are in red. Furthermore, denote by δω := ω(v) − ω the relative angular velocity. The angular part of the Navier-Stokes velocity is locally reconstructed from the velocity field in every lattice, i.e. in every finite element cell T , by means of where x i and t i , for i = 1, 2, 3, 4 are the four nodes and tangential vectors of the lattice and d T = √ 2h T is the diameter of the lattice, see also Figure 3.

Design of the artificial neural network
We will train a deep neural network for determining the coefficients C d , C l , C p and C r . We train two separate neural networks since the input data for C d , C l and C p depends on the angle of attack ψ, while that of C r is invariant to the orientation of the particle. We call these artificial neural networks N and N r . Both take the platelet parameters (L x , L y , L z , α top , α bot ) as input. N further depends on the effective angle of attack ψ P . Alltogether Both neural networks are fully connected feedforward networks with three hidden layers consisting of 50, 20 and 20 neurons in the case of the drag/lift network and 20 neurons each in the case of the rotational torque network. All neurons apart from the output layer are of ReLU type, i.e. using the activation function f (x) = max{x, 0}. Fig. 4 shows the general configuration.

Generation of the training data
Training and test data is obtained by resolved simulations with random sampling of prototypical platelet shapes. Let Ω = {x ∈ R : |x| 2 < 50 µm} \ P be the open ball with radius R = 50 µm around a platelet P . Each platelet P is For training, the Navier-Stokes equations are formulated in the units µm for length, µm · s −1 for the velocity and µg for mass. With the blood viscosity µ = 3 µg · µm −1 · s −1 the Stokes equations are considered. We prescribe zero Dirichlet data on the platelet, v = 0 on ∂P , and set a freestream velocity on the outer boundary ∂Ω \ ∂P . This is either a uniform parallel flow field or a uniform rotational flow field that corresponds to the rotational velocity ω = 2π, both given as Dirichlet data where ψ ∈ [0, 2π] is the relative angle of attack. For v d it holds |v d ψ | = 1 µm · s −1 and in case of the rotational flow it holds |v r ω | = 2πR µm · s −1 , such that it corresponds to a angular velocity of magnitude 2π in counter-clockwise direction around the z-axis.
The training data is generated as follows Algorithm 1 (Generation of the training data) Let N ∈ N be a prescribed number of experiments. For n = 1, 2, . . . , N 1. Generate a random particle P n that satisfies the bounds (7) 2. For four random angles of attack ψ n,i ∈ [0, 2π], i = 1, 2, 3, 4, solve the Stokes equations for the directional Dirichlet data v d ψn,i and compute 1 3. Solve the Stokes equations for the rotational Dirichlet data v d ω and compute the rotational torque Hereby, a set of 4N training data sets (P n , ψ n,i ; F n,i , T n,i ) and N data sets for the rotational configuration (P n ; T n,r ) are generated in an offline phase. Two different networks will be used for these two different settings.
The domain Ω is meshed with hexahedral elements and the finite element discretization is build on equal-order tri-quadratic finite elements for velocity and pressure. The curved boundaries (both the outer boundary and the platelet boundary) are approximated in an isoparametric setup to avoid dominating geometry errors see [29,Sec. 4.2.3]. A very coarse mesh with initially only 12 hexahedras is refined twice around the platelet boundary once globally. The resulting discretization has about 2 000 elements and 60 000 degrees of freedom. Details on the discretization are given in Section 3. The resulting (stationary) discrete finite element formulation is given by For the Stokes equation no transport stabilization must be added.
Given v, p the resulting forces are computed by The units of F and T are Training and test data is computed on an Intel Xeon E5-2640 CPU at 2.40 GHz using 20 parallel threads. A total of 58 500 data sets (46 600 for drag we denote the counter-clockwise rotation of the corresponding vectors by 90 • .  Fig. 5 Visualization of the resolved flow pattern around randomly created particles. The platelets vary in size (Lx × Ly × Lz) and in their convexity, further we vary the angle of attack. The upper row shows three simulations with random particle using different angles of attack. In the lower row, we consider the same particle for each of the three simulations. The two figures on the left correspond to the directional inflow at different angles of attack while the plot on the right corresponds to a simulation with the rotational Dirichlet pattern on the outer boundary of the domain. and lift, 11 900 for measuring the torque) have been generated. The overall computational time for all these 3d simulations was about 9 hours (less than one second for each simulation). All computations are done in Gascoigne 3D [4]. In Figure 5 we show snapshots of three such simulations.

Preparation and normalization of data / training of the neural network
We produce a data set with N entries with random particles. We start by extracting drag, lift and torque according to Algorithm 1. To prepare the input data we encode as much model knowledge as possible. Assume that D, L, T, T r are the vectors containing drag, lift and pitching torque and rotational torque. Then, we define the input data as (component-wise) , l := L 10 sin(2ϕ) , t p := T 2 cos(ϕ) , t r := T r 800 .
These simple relations have been found manually by analyzing the relation of the functional outputs on the different parameters. This scaling reduces the variation of the forces over all experiments to about than 10 − 40% in the case of drag, lift and rotational force. A rescaling of the pitching torque (which has a rather low value) is more difficult since it depends on slight variations of the particle symmetry. The neural network network is implemented in PyTorch [27], using the PyTorch C++ API which has been linked to our finite element framework Gascoigne 3D [4]. The randomly generated data sets originating from the detailed Navier-Stokes simulations are split into 80% serving as training data and the remaining 20% as test data. As loss function we consider square l 2 norm of the error. The training of the two very small networks is accomplished in some few minutes.  Fig. 6 Visualization of the training data. Left: raw values coming from 5 000 experiments (plotted along the x-axis) with randomly generated platelets and random variation of the angle of attack. Right: scaled input data for the neural network according to (18). By prescaling the forces we can reduce the variation to about 20−40%. This remaining dependency of the quantities on the platelet-and flow-parameters will be learned in the artificial neural network.

Testing
To test the accuracy of the trained network we apply it to a set of testing data that was not used in the training of the network. In Figure 7 we show for drag, lift and torque, 250 data points each that have been randomly taken from the test data set such that these data pairs have not been used in training. In the figure, we indicate the exact values as taken from detailed finite element simulations that resolve the particle as large circles and the predicted DNN output as smaller bullets. We observe very good agreement in all three coefficients, best performance in the lift coefficient and highest deviation in the drag.
In Table 1 we indicate the mean (measured in the l 2 -norm) and the maximum error of the network applied to the training data and to the test data. Further, for getting an idea on the generalizability of the approach we also apply the network to additional testing data with random platelets, where at least one of the coefficients (L x , L y , L z , α top , α not ) does not satisfy the bounds specified in (7). We note that such particles are not appearing in the coupled Navier-Stokes particle simulation framework. The average errors appearing in all training and testing data is less than 1%. Maximum relative errors for few single particles reach values up to 4% in the case of the pitching torque, which  Table 1 Accuracy of the neural network model for predicting drag, lift and pitching and rotational torque in percent. We indicate the values for the training data, the test data and the hard test data that consists of data points outside the bounds (7).  Fig. 7 Performance of the neural network in predicting drag, lift and torque coefficients for the flow around randomly created platelets. For 250 random particles each (all have not been used in training the network) we compare the prediction (blue bullets) with the coefficients obtained in a resolved finite element simulation. The coefficients are given in the units of the training reference system described in (17).
is most sensitive with values close to zero. Even if we consider data points that are not within the bounds, average errors are still small although substantial errors are found for single particles.

Application of the neural network
The neural network predicts the coefficients for drag d, lift l, pitching torque t p and rotational torque t r , which are scaled according to (18). While drag, lift and pitching torque depend on the effective angle of attack, the rotational torque is a fixed value that must be predicted only once for each particle. The former three values are recomputed whenever the configuration is changing, i.e. before every advection step.
Algorithm 2 (Neural network / particle / finite element coupling) Let k be the macro time step used to predict the blood flow field, k P := k/M p be the subcycling step for the particle dynamics and N p be the number of particles. For n = 1, 2, . . . iterate 1. Solve the Navier-Stokes equations t n−1 → t n 2. Transfer the local velocities from the finite element mesh to the particle lattice and locally compute the rotational velocity according to (12) 3. For m = 1, . . . , M p subcycle the particle dynamics with step size k P (a) For each particle {P 1 , . . . , P Np }, compute the effective angle of attack accorting to (11) (b) Evaluate the deep neural network for all particles (c) Rescale the coefficients according to (18) and to correct the reference units 2 D = 10 −6 · 70 − 10 cos(2ψ) t L = 10 −6 · 10 sin(2ψ)l T p = 10 −9 · 2 cos(ψ)t p T r = 10 −9 · 800t r (d) Advect all particles and perform collisions according to Remark 1.

Remark 1 (Collisions)
In order to detect and perform collisions we treat particles as spheres with radius L x and use model described in [21]. Since platelets constitutes only small part of the blood volume (less than 1%) collisions between them happen very rarely and this simplification does not affect validity of presented approach. Particles are organized on a lattice mesh. The dimensions of the lattice are generated such that a small number of particles reside in each lattice. This gives a natural way for parallelization and it also helps to keep the communication for performing particle-particle interactions local, compare [24] for details on this appraoch and for a review on further realization techniques.
Step 4 of the algorithm involves the evaluation of the deep neural network. Here, we integrate C++ bindings of the library PyTorch [27] into Gascoigne [4]. All particles are processed at once, such that the evaluation can be performed efficiently in the core of PyTorch. Considering larger networks or a larger number of particles, the use of a CUDA implementation is possible without further effort.
Finally, step 1 of the algorithm requires to solve the Navier-Stokes equations in a finite element framework. The parallel framework that is used in Gascoigne is described in [10,18].

Evaluation of the Navier-Stokes / DNN particle coupling
We study how the different shapes of the particles affect their movement and whether the neural network model is able to give distinguished responses for different particle types, even if the variations of the considered particles are small. In order to do that we examine hydrodynamic forces acting on differently shaped particles. Simulations were performed for five particles (shown in Table 2) representing various shape features (symmetric, asymmetric, convex, convey and combinations).  Table 2 Parameters and shape of 5 test particles. The spatial dimensions are given in µm.
Our domain is a channel of diameter L = 2 mm and infinite length. The schematic geometry of the domain is described in Figure 8. Platelets are variations of ellipsoids with major axes L x × L y × L z with L x ≈ L z ≈ 3 µm and L y ≈ 0.5 µm, for more details see Section 4.1. The inflow data is defined by a time dependent parabolic inflow profile with inflow speed v in = 5 mm/s, namely v in (t) := y All five particles are initially located at y = 0.5176, below the symmetry axis of the velocity profile such that a rotational velocity field attacks the particles. The fluid viscosity is set to µ = 3 mg/mm · s, and particle and fluid density equal ρ = ρ p = 1.06mg/mm 3 . The parameters have been chosen so as to reflect a typical vessel, and realistic blood and platelet properties.
The simulations are carried out with the coupled interaction loop described in Algorithm 2. This means that after each Navier-Stokes step, the fluid velocity is transferred to the particle model and the coupling coefficients drag, lift and torques are updated based on the previously trained neural network. Detailed simulations around different particle shapes only enter the training phase by generating random data sets. Γ in Γ wall particle x y Fig. 8 Spatial configuration of the considered model.

Drag
Angle of incidence Drag coefficient C d particle 1 particle 2 particle 3 particle 4 particle 5 Fig. 9 Drag coefficient as a function of angle of incidence for the five particles defined in Table 2.
Drag is a force acting opposite to the relative motion of the particle moving with respect to a surrounding fluid. Shape-specific drag coefficients present in the literature are usually functions of the particle Reynolds number, angle of incidence and some shape parameters [16,22,44], while drag force itself usually depends on the properties of the fluid and on the size, shape, and speed of the particle.
In Figure 9 the drag coefficient is plotted as a function of the angle of incidence (effective angle of attack) for five considered particles. The first main observation lies in the increase of the drag values when the angle of incidence approaches ψ = π 2 and ψ = 3π 2 so when particle is perpendicular to the flow, which means the biggest cross sectional area with respect to the flow. Correspondingly, the drag decreases when the angle of incidence reaches φ = 0 or φ = 2π so when the cross sectional area gets smaller. Qualitatively, the present results are in good agreement with those issuing from the literature since a similar trend is observed (see e.g. [32]). Furthermore, Figure 9 shows various drag coefficient values for different particles. Particle 2 is characterized by the highest value of the drag. Reasons may be threefold: big size of the particle in comparison to others and hence bigger cross sectional area and higher particle Reynolds number. In contrast, particle 3 is characterized by the lowest value of the drag, which is a result of its small size in comparison to other particles. Particle 1 is an ellipse and serves as a reference. Its drag coefficient is in the middle which is in the line with intuition -particle 1 has intermediate values both in terms of size and convexity/concavity.

Lift
Angle of incidence Lift coefficient C l particle 1 particle 2 particle 3 particle 4 particle 5 Fig. 10 Lift coefficient as a function of angle of incidence for the five particles defined in Table 2.
Lift force on a particle is a result of non-axisymmetric flow field. The pressure distribution on the surface of a particle inclined to the flow direction no longer follows the symmetry of that particle. This gives rise to a lift force due to the displacement of the center of pressure. Lift acts in the direction perpendicular to the fluid velocity and is present when the particles principle axis is inclined to the main flow direction. As in the case of drag, lift coefficient is usually a function of the particle Reynolds number, angle of incidence and some shape parameters [17,25,44], while lift force itself usually depends on the properties of the fluid and on the size, shape, and speed of the particle.
The lift coefficient behaviour at various angles of incidence for five studied particles is presented in Figure 10. The figure shows that the lift coefficient reaches its maximum when the angle of incidence reaches ψ = π 4 or ψ = 5π 4 and its minimum when ψ = 3π 4 or ψ = 7π 4 and is equal to 0 for ψ ∈ {0, π 2 , π, 3π 2 }. These results are consistent with the definition of the lift force and are similar to other studies (see e.g. [26,32]).
Moreover, one can notice that the lift coefficient takes the lowest value for particle 3 and the highest value for particle 2. It results from the difference in surface area, which is small for particle 3 and big for particle 2. Similarly to the drag, the lift of the reference particle 1 is in the middle which also corresponds to the intermediate value of its surface area. Angle of incidence Rotational torque coefficient Cr particle 1 particle 2 particle 3 particle 4 particle 5 Fig. 11 Rotational torque coefficient as a function of angle of incidence for the five particles defined in Table 2. The rotational torque does not depend on the orientation of the particles since it is triggered by the symmetric rotational flow around the particle.
There are two contributions to the rotational motion of the particle. The first is the inherent fluid vorticity, which acts on the particle as a torque due to the resistance on a rotating body. Figure 11 illustrates the rotational torque coefficient plotted as a function of the angle of incidence for five examined particles. One can notice that the magnitude of the coefficients corresponds to the surface area of particle, with particle 2's rotational torque coefficient being the highest, while particle 3 experiencing the smallest rotational torque.
These results are consistent with the definition of the rotational torque and qualitatively are similar to those obtained in the literature (see e.g. [44]).

Pitching torque
Angle of incidence Pitching torque coefficient Cp particle 1 particle 2 particle 3 particle 4 particle 5 Fig. 12 Pitching torque coefficient as a function of angle of incidence for the five particles defined in Table 2.
Since the center of pressure of the total aerodynamic force acting on each particle does not coincide with the particle's center of mass, a pitching torque is generated. This is the second factor that contributes to rotational motion. It accounts for the periodic rotation of the particle around an axis parallel to the flow direction.
In Figure 12 the pitching torque coefficient is plotted as a function of the angle of incidence for all five considered particles. One can notice that the pitching torque coefficient is equal to 0 for ψ = π 2 and ψ = 3π 2 for all five particles, so when particles are perpendicular to the flow. It means that for asymmetric particles, i.e. particle 4 and 5, the pitching torque is 0 when they are set symmetrically with respect to the flow. It may imply that this is their preferred orientation. In case of particle 4 the pitching torque coefficient reaches its minimum when the angle of incidence is ψ = π and its maximum when ψ = 0 or ψ = 2π are reached. It is caused by its asymmetric shape and setting with respect to the direction of the local fluid vorticity, namely particle 4 is convex "at the bottom" and concave "at the top" (for angle of incidence ψ = 0 or ψ = 2π), while the fluid around it is moving clockwise (see Figure 13). For particle 5 the situation is analogous, however it is convex "at the bottom" and concave "at the top". This is consistent with what happens for particle 4 and is reflected on the plot. For the remaining particles 1, 2, 3, the pitching torque is equal or close to 0, which results from their symmetry.  Fig. 13 The pitching torque coefficient Cp depends on particle settings.
In the case of the pitching torque coefficient it is not straightforward to make a comparison between presented trends and those obtained in literature (e.g. [17,25,44]). Most simulations are performed for non-spherical but symmetric particles. Therefore, the discrepancy cannot be easily explained.

Oscillatory translational motion
The translational motion of non-spherical particles is characterized by an oscillatory motion. This is due to the fact that the pressure distribution causes the hydrodynamic forces to work at the center of pressure rather than at the  Table 2.
center of mass. The non-coincidence of the center of pressure and center of mass causes the sustained oscillations (see Figure 14). Moreover, it is observed that every particle is also slowly moving up towards the horizontal axis of symmetry of the domain (all particles start below the axis, see Figure 8). In Figure 15 evolution in time of the aggregated lift of five studied particles is plotted together with the evolution in time of their lift coefficients. One can easily see that the oscillatory motion shown in Figure 14 is a direct consequence of the lift force acting on the particles, while upward motion results from the aggregated lift being positive all the time. The behaviour of the y-velocity is also worth noting (see Figure 16). One can notice that some particles (i.e. 1, 3, 5) decelerate when they are reaching local minimum or maximum. Those local maxima and minima appear for angle of incidence ψ ∈ { π 4 , 3π 4 , 5π 4 , 7π 4 }, so when particles are inclined to the flow direction. Particle 3, the thinnest one, is subjected to the highest deceleration, whereas particles 2 and 4 move more smoothly. Time Lift coefficient and aggregated lift particle 1 particle 2 particle 3 particle 4 particle 5 Fig. 15 Comparison of the lift coefficient for the different particles. In bold lines we show the aggregated lift over time.

Performance of the coupled model for many particles
In this section, we demonstrate the efficiency of the finite element/neural network approach and present numerical results with a multitude of particles. We test the computational effort for the particle model in comparison to the finite element Navier-Stokes discretization. Although some effort has been spent on the multicore implementation based on OpenMP, our implementation is by no means a high performance code. In particular, no use of GPU acceleration within the particle model is applied, neither in the coupling to the neural network model nor in the particle dynamics itself. Both is possible and in parts Time Velocity in y direction particle 1 particle 2 particle 3 particle 4 particle 5 Fig. 16 Velocity in y direction for the five particles defined in Table 2.
already standard in available software packages such like PyTorch C++ [27] or particle dynamics libraries such as LAMMPS [28]. All computations have been carried out on a two-socket system with Intel Xeon E5-2699A v4 processors running at 2.40 Ghz.
We will describe a prototypical blood-flow configuration and discuss the scaling of the implementation with respect to the number of cores. In particular we will investigate the relation between computational effort used in the particle model and in the Navier-Stokes solver. As in Section 5.1 all the parameters have been chosen so as to resemble vessel, blood and platelet properties.

Parallelization
The finite element model is implemented in Gascoigne 3D [4] and outlined in Section 3.1, the discrete systems are approximated with a Newton-Krylow solver using geometric multigrid preconditioning. Basic finite element routines and the linear algebra workflow is partially parallelized using OpenMP, see [10] for details. Since the mesh handling and i/o are not parallelized, substantial speedups are only reached for complex 3d problems.
The particle model is based on a regular lattice mesh that covers the computational domain. The lattice elements of size L h × L h contain the individual particles. Detection of particle-particle collision is limited to those particles that reside in the same lattice element or that belong to directly adjacent elements. This substantially helps to reduce the computational effort which scales quadratically with the number of particles within each lattice element. Hence, we keep L h > 0 small, such that an average of less than 100 particles resides in each element. On the other hand, L h must be chosen large enough to avoid motion of particles accross multiple elements in one time step, i.e. k pd V p ≤ L h , where k pd > 0 is the time step size of the particle model and V p the maximum velocity of the particles. Further, the lattice mesh is basis for parallelization since we can guarantee that no interaction between lattice elements which are separated by a complete layer can take place. We run simulations for the 2D flow in a channel with a local narrowing of 25% which should mimic a stenosed region of a blood channel. Figure 17 displays schematic geometry of the flow domain. The size of the domain, 2 mm × 12 mm, is similar to the dimension of small arteries. The flow is driven by a Dirichlet profile on the inflow boundary Γ in given by

Configuration of the test case
with v in = 5 mm/s. On the wall boundary Γ wall we prescribe no-slip boundary conditions v = 0 and on the outflow boundary Γ out we use the do-nothing outflow condition µ∂ n v − pn = 0, see [14]. The fluid viscosity is set to µ = 3 mg/mm · s. Particle and fluid densities are equal ρ = ρ p = 1.06mg/mm 3 . Due to the fact that platelets constitute less than 1% of the blood volume [19] and the size of the domain we perform simulations with 165 000 particles.
In all numerical examples the temporal step size for solving Navier-Stokes equations is k ns = 0.005 s, while time sub-step for particle advection is k pd = 0.00025 s such that 20 subcycles are computed in each Navier-Stokes step. We update the force coefficients by evaluating the neural network every 10th step (i.e. twice in each Navier-Stokes step). The spatial finite element discretization is based on quadratic elements, with a total of 12819 degrees of freedom.
At the first time step we randomly seed about 165 000 particles distributed over the complete computational domain. Each particle is generated with random properties, i.e. specifying P = (L x , L y , L z , α top , α bot ) by means of the limits indicated in (7) such that the full variety of dimensions and shape is present. Details on the procedure for parametrization of the particles are given in Section 4.1. Then, 10 iterations of the interaction loop shown in Algorithm 2 are performed. Hence, 10 time steps of the Navier-Stokes problem and 200 particle dynamics substeps are performed. Fig. 18 shows the runtime for all 200  Fig. 18 Left: runtime (in seconds) for the coupled Navier-Stokes particle dynamics simulation for an increasing number of cores. Right: parallel speedup for the complete simulation and for the Navier-Stokes finite element simulation and the particle dynamics simulation separately.
iterations. Furthermore we indicate the parallel speedup. These results show that the allocation of computational time to the Navier-Stokes finite element solver and the particle dynamics system is rather balanced. While it is nontrivial to get a reasonable parallel speedup for highly efficient multigrid based finite element simulations (at least for simple 2d problems like this Navier-Stokes testcase), the scaling of the particle dynamics system is superior. These results demonstrate that the number of particles is not the limiting factor for such coupled simulations.
The key feature of our coupled Navier-Stokes particle dynamics scheme is the prediction of the hemodynamical coefficients by means of the previously trained neural network instead of using analytical models, which are not available, or running resolved simulations, which is not feasible for such a large number of particles. In Fig. 19 we give details on the computational time spend in the different parts of the particle dynamics system. Besides advection of the particles, the evaluation of the neural network is dominant, alghouth we update the coefficients in every 10th step only. Here, a more systematic study of the impact of the update frequency should be performed. A further acceleration of the neural network evaluation is possible by using the GPU implementation of PyTorch.
The transition to more realistic 3d problems will substantially increase the effort in both parts, finite elements and particle dynamics. For the Navier-Stokes simulation it has been demonstrated that realistic 3d blood flow situations can be handled in reasonable time, see [8,9,10]. If the number of particles is to be substantially increased, the neural network coupling for estimating hemodynamic coefficients should be realized in a high performance package such as LAMMPS [28] that allows for an efficient GPU implementation.  Fig. 19 Left: runtime (in seconds) for the particle dynamics simulation (all 200 substeps). Right: parallel speedup for the particle advection, handling of particle collisions and the neural network access for predicting the hemodynamical coefficients.

Conclusions
Suspensions of arbitrarily-shaped particles in a fluid are of great importance both in engineering and medical applications. However, the interaction of the non-spherical particles with a fluid flow is a complex phenomenon, even for regularly-shaped particles in the simple fluid flows. The main difficulty lies in determining hydrodynamic forces experienced by a particle due to their strong dependence on both particle shape and its orientation with respect to the fluid flow.
In this paper, a model is successfully derived to simulate the motion of non-spherical particles in a non-uniform flow field, including translation and rotation aspects. The model is designed to reflect platelets in a blood flow, both in terms of particle parameters and fluid configuration. The very good agreement of these results obtained by the coupled finite element / neural network / particle dynamics simulation with state of the art documentations in literature indicates an effectiveness of the presented approach and hence an encouraging potential toward medical applications. Furthermore, the big improvement over usual analytical interaction models is clearly seen as the neural network based model holds for a broad range of different shapes at any orientation. Moreover, using neural network to identify the transmission of forces from fluid to the particles provides a possibility to adopt the model to any desired shape of particle, making this method very promising.
We have further documented details on the scaling of the approach to many particles, which, in 2d blood flow simplifications, matches the typical particle density found for thrombocytes in blood flows. The computational effort is well balanced into the Navier-Stokes finite element part, the particle advection and the evaluation of the neural network.