A finite element / neural network framework for modeling suspensions of non-spherical particles. Concepts and medical applications
NNoname manuscript No. (will be inserted by the editor)
A finite element / neural network framework formodeling suspensions of non-spherical particles
Concepts and medical applications
Martyna Minakowska · Thomas Richter · Sebastian Sager
Received: date / Accepted: date
Abstract
An accurate prediction of the translational and rotational motionof particles suspended in a fluid is only possible if a complete set of cor-relations for the force coefficients of fluid-particle interaction is known. Thepresent study is thus devoted to the derivation and validation of a new frame-work to determine the drag, lift, rotational and pitching torque coefficients fordifferent non-spherical particles in a fluid flow. The motivation for the studyarises from medical applications, where particles may have an arbitrary andcomplex shape. Here, it is usually not possible to derive accurate analyticalmodels for predicting the different hydrodynamic forces. However, consideringfor example the various components of blood, their shape takes an importantrole in controlling several body functions such as control of blood viscosityor coagulation. Therefore, the presented model is designed to be applicableto a broad range of shapes. Another important feature of the suspensions oc-curring in medical and biological applications is the high number of particles.The modelling approach we propose can be efficiently used for simulations ofsolid-liquid suspensions with numerous particles. Based on resolved numericalsimulations of prototypical particles we generate data to train a neural net-work which allows us to quickly estimate the hydrodynamic forces experiencedby a specific particle immersed in a fluid.
Keywords
Non-spherical particles · Finite Element Method · Neural Network
M. Minakowska and T. RichterInstitute of Analysis and Numerics, Otto-von-Guericke University Magdeburg,Universit¨atsplatz 2, 39106 Magdeburg, Germany,E-mail: { martyna.minakowska, thomas.richter } @ovgu.deS. SagerInstitute of Mathematical Optimization, Otto-von-Guericke University Magdeburg,Universit¨atsplatz 2, 39106 Magdeburg, Germany,E-mail: [email protected] a r X i v : . [ phy s i c s . f l u - dyn ] S e p Martyna Minakowska et al.
The prediction of the motion of non-spherical particles suspended in a fluidis crucial for the understanding of natural processes and industrial applica-tions. In such processes, particles can have different shapes and sizes, may bedeformed and can interact with each other. So far, in the majority of scien-tific studies, particulate flow modelling is investigated with the hypothesis ofperfectly spherical particles, thereby eliminating orientation and shape effects.This assumption is very convenient due to its simplicity, the fact that the be-haviour of spheres is well known and the availability of a number of modelsto describe the interaction with fluid flow. The study of suspensions of mul-tiple, irregular-shaped, interacting and deformable particles has received lessattention and still presents a challenge.Particles come in all sort of shapes and sizes, in fact, due to the arbi-trary nature of naturally occurring particles there are an indefinite number ofpossible shapes. On the other hand, there is a common understanding thatparticle shape has a strong influence on the dynamics of NSPS (non-sphericalparticulate systems). These two factors combined makes modelling of NSPS ingeneral way impossible, since for describing the motion of non-spherical parti-cles, detailed information on the fluid dynamic forces acting on such particlesare necessary, but generally not available. Therefore, particular models placeemphasis on different shapes and types of flow. In our work we focus on medi-cal applications, namely on modelling of platelets dynamics under blood flow.Platelets play a main role in the process of blood coagulation and thereforeare of great interest in the modelling of blood flow. The majority of modelscharacterize platelet motion quantitatively and use approaches such as im-mersed boundary method [13], cellular Potts model [39,40,41] and dissipativeparticle dynamics [12,35], treating platelets as points and thus neglecting en-tirely their shape. Some effort has been done to model platelets as rigid two-or three-dimensional spheres or spheroids [23,34,45] but this simplification ofshape 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 par-ticles of irregular shapes, but motivation behind them usually arises fromengineering applications (dispersion of pollutant, pulverized coal combustion,pneumatic transport) [46], where particles are much bigger and usually con-stitute a significant part of the volume of the suspension [33]. But even undergiven very specific circumstances there is no set of correlations of the forcesacting on irregular-shaped particles suspended in a fluid (forces arising fromfluid-particle interaction). Furthermore, in these kinds of models interactionsbetween 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 specificityof the modelled phenomenon. Firstly, platelets, because of their small size(compared to blood vessel), are often modelled quantitatively, ignoring theimportance of shape effects. Secondly, platelets constitute only a very small on-spherical particles in medical flow problems 3 volume of blood what makes most of the engineering applications-driven mod-els, where particles are very dense and often interact with each other, inappro-priate. Thirdly, highly irregular shape of platelets requires new methods forestimating force coefficients in fluid-particle interaction. Fourthly, platelets arenumerous and the evolution of their movement needs to be evaluated quicklyand efficiently.A common approach to the problem of modelling of NSPS consists in de-veloping 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 givinganalytical expressions. Instead, by prototypical simulations we train a neuralnetwork model that takes a couple of parameters describing the particle shape,size and the flow configuration as input and which gives hydrodynamical co-efficients like drag, lift and torque as output. We demonstrate that such amodel 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 requiresthe evaluation of the network for getting updates on these coefficients. Thistwo-step approach with an offline phase for training a network based on theparticle classes under consideration allows for a direct extension to furtherapplications.A similar approach is considered by in [42]. Here, the authors design andtrain a radial basis function network to predict the drag coefficient of non-spherical particles in fluidized beds. Training is based on experimental dataand the network input is the particle’s sphericity and the Reynolds number,covering the Stokes and the intermediate regime. Our approach, based ontraining data generated by detailed simulations of prototypical particles forpredicting drag, lift and torque coefficients could by augmented by includingfurther experimental data.In the following section we will describe prototypical medical applicationswhere such a heterogenous modeling approach can be applied. Then, in Sec-tion 3 we detail the general framework for coupling the Navier-Stokes equa-tion with a discrete particle model. Section 4 introduces the neural networkapproach for estimating the hydrodynamical coefficients and we describe theprocedure for offline training of the network. Then, different numerical testcases are described in Section 5.1 that show the potential of such a heteroge-neous modeling approach. Section 5.2 is devoted to a numerical study com-prising many particles and shows the efficiency of the presented approach. Wesummarize with a short conclusion in Section 6.
Platelets are a vital component of the blood clotting mechanism. They aresmall non-nucleated cell fragments. They have a diameter of approximately2 − µ m, thickness of 0 . µ m, volume of about 7 µ m and a number density of1 . − · µ l − [11] which leads to a volume fraction of only about 10 − : 1. Martyna Minakowska et al.
Still they are a vital component of the blood clotting mechanism. In the reststate platelets shape is discoid, but they have the ability of deforming asa response to various stimuli (chemical and mechanical). They may becomestar shaped (rolling over blood vessels wall to inspect its integrity). Duringthe clotting process they undergo deep morphological changes, from becomingspherical and emitting protuberances (philopodia or pseudopods) which favourmutual aggregation, as well as adhesion to other elements constituting the clotto fully flat, spread stage, to enable wound closure. Thrombocytes constituteapproximately less than 1% of the blood volume, therefore individual plateletshave a negligible effect on blood rheology [19].Due to a significant effect of the particles shape on their motion and theirpractical importance in industrial applications, the non-sphericity started at-tracting attention in the modelling and simulations of particles transport influid flows [17,22,44]. Unfortunately, it is not possible to consider each shapein the implementation of numerical methods because of non-existence of asingle approach describing accurately the sizes and shapes of non-sphericalparticles. Spheres can be described by a single characteristic value, i.e. thediameter, whereas non-spherical particles require more parameters. Even veryregular shapes need at least two parameters. Moreover, the particles may havevarying orientation with respect to the flow, what makes the description oftheir behaviour even more difficult. Even though several methods for shapeparametrization and measurement have been suggested, none has won greateracceptance. One of most commonly used shape factor is a sphericity which wasfirstly introduced in [37] and defined as the ratio between the surface area of asphere with equal volume as the particle and the surface area of the consideredparticle. Then the drag coefficient for non-spherical particle is estimated byusing correlations for spherical particles and modified to take into account thesphericity factor [43]. Models using sphericity as a shape factor give promis-ing results when restricted to non-spherical particles with aspect ratios β lessthan 1.7 [7], where β = L/D with L and D being length and diameter ofthe considered particle. For particles having extreme shapes and those havinglittle resemblance to a sphere, the sphericity concept fails to produce satis-fying quantitative results [6]. In general, the lower the sphericity, the pooreris the prediction. Also, the same value of the sphericity may be obtained fortwo very varying shapes whose behaviour in the flow is different. Moreover, thesphericity does not take the orientation into account. In order to introduce ori-entation dependency in drag correlations, some researchers use two additionalfactors: the crosswise sphericity and lengthwise sphericity [16]. Most of thesecorrelations employ also dependency on particle Reynolds number defined as Re p = ρ ¯ ud eq /µ, where ρ and µ are fluid density and viscosity, ¯ u = u f − u p is the velocityof the particle relative to the fluid velocity and d eq is the equivalent particlediameter, i.e. the diameter of a sphere with the same volume as the consideredparticle in order to include the importance of fluid properties. on-spherical particles in medical flow problems 5 The shape factor concept may be described as an attempt to define asingle correlation for drag for all shapes and orientations. Another approachappeared as an alternative consisting in obtaining drag coefficient expressionsfor a fixed shape and any orientations: the drag coefficient is determined attwo extreme angles of incidence (0 ◦ and 90 ◦ ) from existing correlations whichare then linked by some functions resulting in the whole range of angles ofincidence for non-spherical particles [30]. However, besides drag force, non-spherical particles are associated with orientation and shape induced lift alongwith pitching and rotational torques. H¨olzer and Sommerfeld [17] investigateda few different shapes of non-spherical particles at different flow incident an-gles using the Lattice Boltzmann method to simulate the flow around theparticle. 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 bya direct numerical simulation (DNS) carried out with an immersed boundarymethod. Those correlations employ particle Reynolds number, angle of inci-dence and some shape-related coefficients. Ouchene et al. [25] determined forcecoefficients depending on particle Reynolds number, aspect ratio and angle ofincidence by fitting the results extracted from DNS computations of the flowaround prolate ellipsoidal particles. Discrete element methods (DEM) coupledwith computational fluid dynamics (CFD) has been recognized as a promisingmethod to meet the challenges of modelling of NSPS [47,48]. DEM is a numer-ical approach for modelling a large number of particles interacting with eachother. The simplest computational sequence for the DEM typically proceedsby solving the equations of motion, while updating contact force histories as aconsequence of contacts between different discrete elements and/or resultingfrom contacts with model boundaries. It is designed to deal with very densesuspensions, where contacts between particles are very common and play a keyrole in determining the motion of particles, see [46] for an extensive overviewof DEM.Obviously, for a particle with a certain specific shape, the general ex-pressions derived from the first factor shape approach tend to be less ac-curate than the specialized one for that shape, but the efficiency of interpola-tions/extrapolations to the various shapes to provide the general expression isan attractive perspective on engineering applications. On the other hand, par-ticles occurring in biological processes are usually very numerous. Therefore,an effective method not only has to be accurate but also efficient in terms ofcomputational time.To overcome the aforementioned limitations in modelling of NSPS we em-ploy the recently trending approach and use machine learning to design amethod that enables us to model the behaviour of suspensions of particlesof an arbitrary shape while maintaining at the same time the accuracy ofshape-specific models. We also place an emphasis on the computational effi-ciency as usually there are plenty of particles involved in medical processesand engineering problems. Martyna Minakowska et al.
This section describes the general numerical framework for suspensions of par-ticles in a Navier-Stokes fluid. The discretization of the Navier-Stokes equa-tions is realized in the finite element toolbox
Gascoigne 3D and outlined inSection 3.1. Then, in Section 3.2 we describe a very simple model for themotion of the particles.3.1 Fluid dynamicsConsider a finite time interval I = [0 , T ] and a bounded domain Ω ∈ R d for d ∈ { , } . We assume incompressibility of fluid, which is modelled by theNavier-Stokes equations that take the form ρ (cid:0) ∂ t v + ( v · ∇ ) v (cid:1) − div σ ( v , p ) = 0 , div v = 0 , (1)where v denotes the fluid velocity, σ is the Cauchy stress tensor σ ( v , p ) = ρν ( ∇ v + ∇ v T ) − p I ,p the pressure, ρ the fluid mass density and ν the kinematic viscosity. Fluiddensity and viscosity are assumed to be nonnegative and constant.The fluid boundary is split into an inflow boundary Γ in , an outflow bound-ary Γ out and rigid no-slip wall boundaries Γ wall . On the inflow and walls weimpose Dirichlet boundary conditions while on the outflow we apply the do-nothing condition (see e.g. [15]) v = v in on Γ in × I, v = 0 on Γ wall × I, ( ρν ∇ v − p I ) n = 0 on Γ out × I, (2)where v in is prescribed inflow-profile and n is the outward unit normal vector. Discretization
For temporal discretization of the Navier-Stokes equations weintroduce a uniform partitioning of the interval I = [0 , T ] into discrete steps0 = t < t < · · · < t N = T, k := t n − t n − . By v n := v ( t n ) and p n := p ( t n ) we denote the approximations at time t n . Weuse a shifted version of the Crank-Nicolson time discretization scheme whichis second order accurate and which has preferable smoothing properties ascompared to the standard version, see [14, Remark 1], i.e. ρ (cid:18) v n − v n − k + θ ( v n · ∇ ) v n + (1 − θ )( v n − · ∇ ) v n − (cid:19) − θ div σ ( v n , p n ) − (1 − θ ) div σ ( v n , p n ) = 0 , (3) on-spherical particles in medical flow problems 7 where, typically, θ = k . For spatial discretization we denote by Ω h a quadri-lateral (or hexahedral) finite element mesh of the domain Ω that satisfiesthe usual regularity assumptions required for robust interpolation estimates,see [29, Section 4.2]. Adaptive meshes are realized by introducing at most onehanging node per edge. Discretization is based on second order finite elementsfor pressure and velocity. To cope with the lacking inf-sup stability of this equalorder finite element pair we stabilize with the local projection method [2]. Localprojection terms are also added to stabilize dominating transport [3]. Finally,velocity v n ∈ [ V h ] and pressure p n ∈ V h (where we denote by V h the space ofbi-quadratic finite elements on the quadrilateral mesh) are given as solutionto (cid:0) ρ v n , φ h (cid:1) Ω + kθ (cid:0) ρ ( v n · ∇ ) v n , φ h (cid:1) Ω + kθ (cid:0) ρν ( ∇ v n + ∇ v Tn ) , ∇ φ (cid:1) Ω − k (cid:0) p n , div φ (cid:1) Ω + k (cid:0) ρ div v n , ξ h (cid:1) Ω + k (cid:88) K ∈ Ω h α K (cid:0) ∇ ( p n − π h p n ) , ∇ ( ξ h − π h ξ h ) (cid:1) K + k (cid:88) K ∈ Ω h α K (cid:0) ( v n · ∇ )( v n − π h v n ) , ( v n · ∇ )( φ h − φ ξh ) (cid:1) K = (cid:0) ρ v n − , φ h (cid:1) Ω − k (1 − θ ) (cid:0) ρ ( v n − · ∇ ) v n − , φ h (cid:1) Ω − k (1 − θ ) (cid:0) ρν ( ∇ v n − + ∇ v Tn − ) , ∇ φ (cid:1) Ω ∀ ( φ h , ξ h ) ∈ [ V h ] × V h , (4)where the stabilization parameters are element-wise chosen as [5] α K = α (cid:18) νh K + (cid:107) v (cid:107) L ∞ ( K ) h K + 1 k (cid:19) − , where h K = diam( K ) is the diameter of the element K . Usually we choose α = 0 .
1. By π h : V h → V (1) h we denote the interpolation into the space ofbi-linear elements on the same mesh Ω h . Solution of the discretized problem
Discretization by means of (4) gives riseto a large system of nonlinear algebraic equation which we approximate bya Newton scheme based on the analytical Jacobian of (4). The resulting lin-ear systems are solved by a GMRES iteration (Generalized minimal residualmethod [31]), preconditioned by a geometric multigrid solver [1]. As smootherwe employ a Vanka-type iteration based on the inversion of the submatricesbelonging to each finite element cell. These local 27 ×
27 (108 ×
108 in 3d) ma-trices are inverted exactly. Essential parts of the complete solution frameworkare parallelized using OpenMP, see [10].3.2 Particle dynamicsThe particles suspended in the fluid are described as rigid bodies and theirdynamics is driven by the hydrodynamical forces of the flow. Each particle P Martyna Minakowska et al. with center of mass x P , velocity V P and angular velocity Ω P is governed byNewton’s law of motion m P X (cid:48)(cid:48) p ( t ) = F ( v , p ; P ) := − (cid:90) ∂P σ ( v , p ) n P d s J P Ω (cid:48) P ( t ) = T ( v , p ; P ) := (cid:90) ∂P ( x − x P ) × (cid:0) σ ( v , p ) n P (cid:1) d s, (5)where m P is the particle’s mass, and J P its moment of inertia given by J P = diag (cid:110) ρ P (cid:90) P (cid:16) | x − x P | − ( x i − [ x P ] i ) (cid:17) d x ; i = 1 , , (cid:111) , with the (uniform) particle density ρ P . n P is the unit normal vector on theparticle boundary facing into the fluid.A resolved simulation is out of bounds due to the large number of plateletsand in particular due to the discrepancy in particle diameter (about 10 − m)versus vessel diameter (about 10 − m). Instead, we consider all platelets to bepoint-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 onthe shape and the size of the particles but also on their relative orientation andmotion in the velocity field of the blood. Since the relative velocities (bloodvs. particles) are very small the interaction lies within the Stokes regime witha linear scaling in terms of the velocity. The deep neural network will predictcoefficients 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 F P = C d ( P, ψ P )( v − V p ) + C l ( P, ψ P )( v − V p ) ⊥ T P = C p ( P, ψ P ) | v − V p | + C r ( P, ψ P )( ω − Ω p ) , (6)where P = ( L x , L y , L z , α top , α bot ) describes the particle shape and where ψ P isthe relative angle of attack which depends on the particle orientation but alsoon 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 basedon detailed numerical simulations using random particles in random configu-rations. By ( v − V p ) ⊥ we denote the flow vector in lift-direction, orthogonal tothe main flow direction. In 3d configurations, two such lift coefficients must betrained. Here we will however only consider 2d simplifications with one dragand one lift direction. on-spherical particles in medical flow problems 9 In this section we describe the neural network model for coupling the Navier-Stokes equation with a suspension of non-spherical particles. The differenthydrodynamical coefficients will be taken from a neural network, which istrained in an offline phase. Training data is achieved by resolved Navier-Stokessimulations using prototypical particles with random parameters.The setting investigated in this work carries several special characteristicsthat differ from industrial applications. – The particle density is very small - about 1 . − . · gl − and the par-ticle and fluid densities are similar (the average density of whole blood fora human is about 1 . · gl − ). Blood contains about 200 000 – 400 000platelets per mm summing up to less than 1% of the overall blood vol-ume [38]. Hence we neglect all effects of the particles onto the fluid. Thissimplification is possible since we only model the platelets as rigid parti-cles. Effects of the red blood cells, much larger and appearing in greaterquantity, can be integrated by means of a non-Newtonian rheology. – The particle Reynolds numbers are very small (with order of magnitudeabout 10 − or less) such that we are locally in the Stokes regime. This ismainly due to the smallness of the platelets (diameter approximately 3 µ m)and the small flow velocities at (bulk) Reynolds numbers ranging from 50to 1 000 depending on the specific vessel under investigation. We focus oncoronary vessels with a diameter around 2 mm and with Reynolds numberabout 200. – The platelets have a strongly non-spherical, disc-like shape. Their shapeand size underlies a natural variation. Furthermore, under activation, theparticles will take a spherical shape.Instead of deriving analytical models for the transmission of forces from thefluid to the particles, we develop a neural network for the identification ofdrag, lift and torque coefficients based on several parameters describing theshape and the size of the platelets and the individual flow configuration.4.1 Parametrization of the plateletsWe model the platelets as variations of an ellipsoid with major axes L x × L y × L z with L x ≈ L z ≈ µ m and L y ≈ . µ m. In y -direction upper α top and lower α bot semi-ellipsoids are modified to give them a more or less concave or convexshape. 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 ofthe levelset function Φ (cid:0) P ; x, y, z (cid:1) = 1 − R ( P ) − (cid:16) α + (1 − α ) R ( P ) (cid:17) − (cid:18) yL y (cid:19) , − − . . − − . . α top = 0 . , α bot = 0 . α top = 1 . , α bot = 1 . − − . . − − . . α top = 1 . , α bot = 0 . α top = 0 . , α bot = 1 . Fig. 1
Prototypical templates for variations of the parameters L x , L y , L z (length) and α top , α bot (shape). Within the blood flow, particles can appear at all angles of attack φ . where we define R ( P ) := (cid:18) xL x (cid:19) + (cid:18) zL z (cid:19) and α := (cid:40) α top y > = 0 α bot y < . We assume that all parameters L x , L y , L z , α top , α bot are normally distributedwith means indicated above and with standard derivation 0 . L x , L y , L z and 0 . α top , α bot . We drop particles thatexceed the bounds L x , L z ∈ [2 . µ m , . µ m] , L y ∈ [0 . µ m , µ m] , α top , α bot ∈ [0 . , . (7)In Fig. 1 we show some typical shapes of the platelets.Next, we indicate mass, center of mass and moment of inertia for a parametrizedparticle P = ( L x , L y , L z , α top , α bot ). Unless otherwise specified, all quantitiesare given in µm and g . The mass of a particle P is approximated by m ( P ) := ρ P L x L y L z . (cid:0) √ α bot + 0 . (cid:112) α top + 0 . (cid:1) . (8)This approximation is based on a weighted one-point Gaussian quadraturerule. It is accurate with an error of at most 2% for all α ∈ [0 . , P is given by m x ( P ) = m z ( P ) = 0 , m y ( P ) = L x L y L z π m ( P ) (cid:0) α top − α bot (cid:1) . (9) on-spherical particles in medical flow problems 11 \ 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 ori-entation. Right: the black arrow δ v = v − V is the velocity acting on the platelet. By ψ := ∠ ( δ v , e x ) − ϕ we denote the effective angle of attack. The moment of inertia in the x/y plane (the only axis of rotation that wewill consider in the 2d simplification within this work) is given by I z ( P ) = ρ P (cid:90) P (cid:0) x + y ) d( xyz ) ≈ L x L y L z ρ P π (cid:16) .
24 + 0 . α top + α bot ) + 0 . α top + α bot ) (cid:17) , (10)which is accurate up to an error of at most 1% for all α ∈ [0 . ,
2d Simplification
To start with, we apply a two dimensional simplification ofthe problem by assuming that the blood vessel is a layer of infinite depth (in z -direction) and that it holds v = 0 for the blood velocity and V = 0 for allparticles. 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, rotationis restricted to rotation around the x/y -plane. Hence, Ω = (0 , , ω ) is describedby a scalar component. A complete particle is then described by P = ( L x , L y , L z , α top , α bot ; X , ϕ, V , Ω ) , 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 Ω theangular velocity w.r.t. the z -axis rotation.To describe the forces acting on the particle suspended in the Navier-Stokesfluid we denote by δ v := v − V the effective velocity vector, i.e. the relativevelocity that is acting on the platelets. By ψ := ∠ ( e x , v − V ) − ϕ we denotethe 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 iscomputed as ψ := ∠ ( e x , v − V ) − ϕ. (11) v ( x i ) the Navier-Stokes veloc-ity at node x i , by v ( t i ) the unit-tangentialvector at this node in counter-clockwise ori-entation. With (cid:104) v i , t i (cid:105) we denote the velocitycontribution in tangential direction. Positivevalues are indicated in orange (at node x ),negative contributions are in red. Fig. 3
Approximation of the rotational velocity on an element T . Furthermore, denote by δω := ω ( v ) − ω the relative angular velocity. Theangular part of the Navier-Stokes velocity is locally reconstructed from thevelocity field in every lattice, i.e. in every finite element cell T , by means of ω T = 12 d T (cid:88) i =1 (cid:104) v ( x i ) , t i (cid:105) , (12)where x i and t i , for i = 1 , , , d T = √ h T is the diameter of the lattice, see also Figure 3.4.2 Design of the artificial neural networkWe 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 invariantto 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 N : R → R , N r : R → R . Both neural networks are fully connected feedforward networks with threehidden layers consisting of 50, 20 and 20 neurons in the case of the drag/liftnetwork and 20 neurons each in the case of the rotational torque network. Allneurons apart from the output layer are of ReLU type, i.e. using the activationfunction f ( x ) = max { x, } . Fig. 4 shows the general configuration.4.3 Generation of the training dataTraining and test data is obtained by resolved simulations with random sam-pling of prototypical platelet shapes. Let Ω = { x ∈ R : | x | < µ m } \ P bethe open ball with radius R = 50 µ m around a platelet P . Each platelet P is on-spherical particles in medical flow problems 13 Fig. 4
Architecture for the neural networks N (left) for predicting drag, lift and pitchingtorque and N r (right) for predicting the rotational torque. Both networks are fully connectedfeedforward networks with 3 hidden layers, having 50/20/20 and 20/20/20 neurons each. constructed by taking normally distributed values for ( L x , L y , L z , α top , α bot )as indicated in Section 4.1. Since only relative velocity and relative orientationmatter, the platelets are fixed in the origin with V = 0 and Ω = 0 at angle ϕ = 0.For training, the Navier-Stokes equations are formulated in the units µ mfor length, µ m · s − for the velocity and µ g for mass. With the blood viscosity µ = 3 µ g · µ m − · s − the Stokes equations − µ∆ v + ∇ p = 0 , div v = 0 , (13)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 auniform parallel flow field or a uniform rotational flow field that correspondsto the rotational velocity ω = 2 π , both given as Dirichlet data v = 0 on ∂P, v dψ := cos( ψ )sin( ψ )0 or v rω := 2 π − yx on ∂Ω \ ∂P, (14)where ψ ∈ [0 , π ] is the relative angle of attack. For v d it holds | v dψ | =1 µ m · s − and in case of the rotational flow it holds | v rω | = 2 πR µ m · s − , suchthat it corresponds to a angular velocity of magnitude 2 π in counter-clockwisedirection around the z -axis.The training data is generated as follows Algorithm 1 (Generation of the training data)
Let N ∈ N be a pre-scribed number of experiments. For n = 1 , , . . . , N
1. Generate a random particle P n that satisfies the bounds (7)
2. For four random angles of attack ψ n,i ∈ [0 , π ] , i = 1 , , , , solve theStokes equations for the directional Dirichlet data v dψ n,i and compute D n,i := (cid:90) ∂P (cid:16) µ ∇ v − pI (cid:17) n · v dψ n,i d s L n,i := (cid:90) ∂P (cid:16) µ ∇ v − pI (cid:17) n · v d, ⊥ ψ n,i d s T n,i := (cid:90) ∂P ( x − x P ) ⊥ · (cid:16) µ ∇ v − pI (cid:17) n d s
3. Solve the Stokes equations for the rotational Dirichlet data v dω and computethe rotational torque T n,r := (cid:90) ∂P ( x − x P ) ⊥ · (cid:16) µ ∇ v − pI (cid:17) n d s. Hereby, a set of 4 N training data sets ( P n , ψ n,i ; F n,i , T n,i ) and N data setsfor 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 elementdiscretization is build on equal-order tri-quadratic finite elements for veloc-ity and pressure. The curved boundaries (both the outer boundary and theplatelet boundary) are approximated in an isoparametric setup to avoid dom-inating geometry errors see [29, Sec. 4.2.3]. A very coarse mesh with initiallyonly 12 hexahedras is refined twice around the platelet boundary once glob-ally. The resulting discretization has about 2 000 elements and 60 000 degreesof freedom. Details on the discretization are given in Section 3. The resulting(stationary) discrete finite element formulation is given by ρν (cid:0) ∇ v , ∇ φ (cid:1) − (cid:0) p, div φ (cid:1) = 0 ∀ φ ∈ [ V h ] (cid:0) div v , ξ (cid:1) + (cid:88) T ∈ Ω h h T ν (cid:0) ∇ ( p − π h p ) , ∇ ( ξ − π h ξ ) (cid:1) T = 0 ∀ ξ ∈ V h . (15)For the Stokes equation no transport stabilization must be added.Given v , p the resulting forces are computed by F = (cid:90) ∂P (cid:0) µ ∇ v − pI (cid:1) n d o, T = (cid:90) ∂P ( x − x P ) ⊥ · (cid:0) µ ∇ v − pI (cid:1) n d o. (16)The units of F and T are[ F ] = µ m · µ gs , [ T ] = µ m · µ gs . (17)Training and test data is computed on an Intel Xeon E5-2640 CPU at 2.40GHz using 20 parallel threads. A total of 58 500 data sets (46 600 for drag By ( x − x P ) ⊥ and v d, ⊥ ψ n,i we denote the counter-clockwise rotation of the correspondingvectors by 90 ◦ .on-spherical particles in medical flow problems 15 Fig. 5
Visualization of the resolved flow pattern around randomly created particles. Theplatelets vary in size ( L x × L y × L z ) and in their convexity, further we vary the angle ofattack. The upper row shows three simulations with random particle using different anglesof 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 attackwhile the plot on the right corresponds to a simulation with the rotational Dirichlet patternon the outer boundary of the domain. and lift, 11 900 for measuring the torque) have been generated. The overallcomputational time for all these 3d simulations was about 9 hours (less thanone second for each simulation). All computations are done in Gascoigne 3D [4].In Figure 5 we show snapshots of three such simulations. We produce a data set with N entries with random particles. We start byextracting drag, lift and torque according to Algorithm 1. To prepare the inputdata 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) d := D −
10 cos(2 ϕ ) , l := L
10 sin(2 ϕ ) , t p := T ϕ ) , t r := T r . (18)These simple relations have been found manually by analyzing the relation ofthe functional outputs on the different parameters. This scaling reduces thevariation of the forces over all experiments to about than 10 −
40% in the caseof drag, lift and rotational force. A rescaling of the pitching torque (which hasa rather low value) is more difficult since it depends on slight variations of theparticle symmetry.The neural network network is implemented in
PyTorch [27], using the
PyTorch C++ API which has been linked to our finite element frameworkGascoigne 3D [4]. The randomly generated data sets originating from thedetailed Navier-Stokes simulations are split into 80% serving as training dataand the remaining 20% as test data. As loss function we consider square l norm of the error. The training of the two very small networks is accomplishedin some few minutes. . − . − . . . . . . . . −
10 1 . . . . t o r q u e p i t c h i n g t o r q u e li f t d r ag 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 theangle of attack. Right: scaled input data for the neural network according to (18). By pre-scaling the forces we can reduce the variation to about 20 − To test the accuracy of the trained network we apply it to a set of testing datathat 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 thetest data set such that these data pairs have not been used in training. In thefigure, we indicate the exact values as taken from detailed finite element simu-lations that resolve the particle as large circles and the predicted DNN outputas 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 -norm) and the max-imum 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 alsoapply the network to additional testing data with random platelets, where atleast one of the coefficients ( L x , L y , L z , α top , α not ) does not satisfy the boundsspecified in (7). We note that such particles are not appearing in the coupledNavier-Stokes particle simulation framework. The average errors appearing inall training and testing data is less than 1%. Maximum relative errors for fewsingle particles reach values up to 4% in the case of the pitching torque, which on-spherical particles in medical flow problems 17training data test data generalizationavg max avg max avg maxDrag 0 .
40% 1 .
44% 0 .
40% 1 .
43% 1 .
35% 6 . .
34% 1 .
32% 0 .
34% 1 .
33% 2 .
22% 19 . .
93% 3 .
75% 0 .
95% 3 .
87% 6 .
62% 55 . .
45% 1 .
57% 0 .
45% 1 .
50% 2 .
21% 11 . Table 1
Accuracy of the neural network model for predicting drag, lift and pitching androtational torque in percent. We indicate the values for the training data, the test data andthe hard test data that consists of data points outside the bounds (7). rotational torque 2502001501005001000800600400 pitching torque 25020015010050010 − − − Fig. 7
Performance of the neural network in predicting drag, lift and torque coefficientsfor the flow around randomly created platelets. For 250 random particles each (all havenot been used in training the network) we compare the prediction (blue bullets) with thecoefficients obtained in a resolved finite element simulation. The coefficients are given in theunits of the training reference system described in (17). is most sensitive with values close to zero. Even if we consider data points thatare not within the bounds, average errors are still small although substantialerrors are found for single particles.4.4 Application of the neural networkThe 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 rotationaltorque is a fixed value that must be predicted only once for each particle. Theformer 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 thesubcycling step for the particle dynamics and N p be the number of particles.For n = 1 , , . . . iterate1. Solve the Navier-Stokes equations t n − (cid:55)→ t n
2. Transfer the local velocities from the finite element mesh to the particlelattice 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 , . . . , P N p } , compute the effective angle of attackaccorting to (11)(b) Evaluate the deep neural network for all particles ( L x , L y , L z , α top , α bot , ψ ) i (cid:55)→ ( d , l , t p , t r ) i i = 1 , . . . , N p (c) Rescale the coefficients according to (18) and to correct the referenceunits D = 10 − · (cid:0) −
10 cos(2 ψ ) (cid:1) t L = 10 − ·
10 sin(2 ψ ) lT p = 10 − · ψ ) t p T r = 10 − · t r (d) Advect all particles and perform collisions according to Remark 1.Remark 1 (Collisions) In order to detect and perform collisions we treat par-ticles as spheres with radius L x and use model described in [21]. Since plateletsconstitutes only small part of the blood volume (less than 1%) collisions be-tween them happen very rarely and this simplification does not affect validityof presented approach.Usually we choose M p = 100 subcycling iterations within each macro step. Forfurther acceleration, steps of the inner loop can be skipped in mostof these inner iterations and it will be sufficient to recalculate the couplingcoefficients in approximately every tenth step. Remark 2 (Parallelization)
Steps are parallelized using OpenMP.Particles are organized on a lattice mesh. The dimensions of the lattice aregenerated such that a small number of particles reside in each lattice. Thisgives a natural way for parallelization and it also helps to keep the commu-nication for performing particle-particle interactions local, compare [24] fordetails 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 performedefficiently in the core of PyTorch. Considering larger networks or a largernumber of particles, the use of a CUDA implementation is possible withoutfurther effort.Finally, step 1 of the algorithm requires to solve the Navier-Stokes equa-tions in a finite element framework. The parallel framework that is used inGascoigne is described in [10,18]. While the blood flow configuration is based in mm and g, the particle configuration isset up in µ m and µ g.on-spherical particles in medical flow problems 19 Particle Shape Lx Ly Lz α α 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. Theschematic geometry of the domain is described in Figure 8. Platelets are vari-ations of ellipsoids with major axes L x × L y × L z with L x ≈ L z ≈ µ m and L y ≈ . µ m, for more details see Section 4.1. The inflow data is defined bya time dependent parabolic inflow profile with inflow speed v in = 5 mm / s,namely v in ( t ) := y (2 − y )2 t · v in . All five particles are initially located at y = 0 . µ = 3 mg / mm · s, and particle and fluid densityequal ρ = ρ p = 1 . / mm . The parameters have been chosen so as to reflecta typical vessel, and realistic blood and platelet properties.The simulations are carried out with the coupled interaction loop describedin Algorithm 2. This means that after each Navier-Stokes step, the fluid ve-locity 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 trainingphase by generating random data sets. Γ in Γ wall particle xy Fig. 8
Spatial configuration of the considered model. π π π π . · − Angle of incidence D r ag c o e ffi c i e n t C d particle 1particle 2particle 3particle 4particle 5 Fig. 9
Drag coefficient as a function of angle of incidence for the five particles defined inTable 2.
Drag is a force acting opposite to the relative motion of the particle movingwith respect to a surrounding fluid. Shape-specific drag coefficients present inthe literature are usually functions of the particle Reynolds number, angle ofincidence and some shape parameters [16,22,44], while drag force itself usuallydepends on the properties of the fluid and on the size, shape, and speed of theparticle.In Figure 9 the drag coefficient is plotted as a function of the angle ofincidence (effective angle of attack) for five considered particles. The first mainobservation lies in the increase of the drag values when the angle of incidence on-spherical particles in medical flow problems 21 approaches ψ = π and ψ = π so when particle is perpendicular to theflow, which means the biggest cross sectional area with respect to the flow.Correspondingly, the drag decreases when the angle of incidence reaches φ = 0or φ = 2 π so when the cross sectional area gets smaller. Qualitatively, thepresent results are in good agreement with those issuing from the literaturesince a similar trend is observed (see e.g. [32]).Furthermore, Figure 9 shows various drag coefficient values for differentparticles. Particle 2 is characterized by the highest value of the drag. Reasonsmay be threefold: big size of the particle in comparison to others and hencebigger cross sectional area and higher particle Reynolds number. In contrast,particle 3 is characterized by the lowest value of the drag, which is a resultof its small size in comparison to other particles. Particle 1 is an ellipse andserves as a reference. Its drag coefficient is in the middle which is in the linewith intuition - particle 1 has intermediate values both in terms of size andconvexity/concavity. π π π π − − · − Angle of incidence L i f t c o e ffi c i e n t C l particle 1particle 2particle 3particle 4particle 5 Fig. 10
Lift coefficient as a function of angle of incidence for the five particles defined inTable 2.
Lift force on a particle is a result of non-axisymmetric flow field. Thepressure distribution on the surface of a particle inclined to the flow directionno longer follows the symmetry of that particle. This gives rise to a lift forcedue to the displacement of the center of pressure. Lift acts in the directionperpendicular to the fluid velocity and is present when the particles principleaxis 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 andsome shape parameters [17,25,44], while lift force itself usually depends onthe 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 studiedparticles is presented in Figure 10. The figure shows that the lift coefficientreaches its maximum when the angle of incidence reaches ψ = π or ψ = π andits minimum when ψ = π or ψ = π and is equal to 0 for ψ ∈ { , π , π, π } .These results are consistent with the definition of the lift force and are similarto other studies (see e.g. [26,32]).Moreover, one can notice that the lift coefficient takes the lowest value forparticle 3 and the highest value for particle 2. It results from the differencein surface area, which is small for particle 3 and big for particle 2. Similarlyto the drag, the lift of the reference particle 1 is in the middle which alsocorresponds to the intermediate value of its surface area. π π π π . . . . · − Angle of incidence R o t a t i o n a l t o r q u ec o e ffi c i e n t C r particle 1particle 2particle 3particle 4particle 5 Fig. 11
Rotational torque coefficient as a function of angle of incidence for the five particlesdefined in Table 2. The rotational torque does not depend on the orientation of the particlessince it is triggered by the symmetric rotational flow around the particle.
There are two contributions to the rotational motion of the particle. Thefirst is the inherent fluid vorticity, which acts on the particle as a torque dueto the resistance on a rotating body.Figure 11 illustrates the rotational torque coefficient plotted as a functionof the angle of incidence for five examined particles. One can notice that themagnitude of the coefficients corresponds to the surface area of particle, with on-spherical particles in medical flow problems 23 particle 2’s rotational torque coefficient being the highest, while particle 3experiencing the smallest rotational torque.These results are consistent with the definition of the rotational torque andqualitatively are similar to those obtained in the literature (see e.g. [44]). π π π π − − . . · − Angle of incidence P i t c h i n g t o r q u ec o e ffi c i e n t C p particle 1particle 2particle 3particle 4particle 5 Fig. 12
Pitching torque coefficient as a function of angle of incidence for the five particlesdefined in Table 2.
Since the center of pressure of the total aerodynamic force acting on eachparticle does not coincide with the particle’s center of mass, a pitching torqueis generated. This is the second factor that contributes to rotational motion.It accounts for the periodic rotation of the particle around an axis parallel tothe flow direction.In Figure 12 the pitching torque coefficient is plotted as a function ofthe angle of incidence for all five considered particles. One can notice thatthe pitching torque coefficient is equal to 0 for ψ = π and ψ = π for allfive particles, so when particles are perpendicular to the flow. It means thatfor asymmetric particles, i.e. particle 4 and 5, the pitching torque is 0 whenthey are set symmetrically with respect to the flow. It may imply that this istheir preferred orientation. In case of particle 4 the pitching torque coefficientreaches its minimum when the angle of incidence is ψ = π and its maximumwhen ψ = 0 or ψ = 2 π are reached. It is caused by its asymmetric shapeand setting with respect to the direction of the local fluid vorticity, namelyparticle 4 is convex “at the bottom” and concave “at the top” (for angle ofincidence ψ = 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 “atthe bottom” and concave “at the top”. This is consistent with what happensfor 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.
A positive C p B negative C p Fig. 13
The pitching torque coefficient C p depends on particle settings. In the case of the pitching torque coefficient it is not straightforward tomake a comparison between presented trends and those obtained in litera-ture (e.g. [17,25,44]). Most simulations are performed for non-spherical butsymmetric particles. Therefore, the discrepancy cannot be easily explained.
The translational motion of non-spherical particles is characterized by an os-cillatory motion. This is due to the fact that the pressure distribution causesthe hydrodynamic forces to work at the center of pressure rather than at the . . . . . . . . . . Time P o s i t i o n i n y d i r ec t i o n particle 1particle 2particle 3particle 4particle 5 Fig. 14
Position in y direction for the five particles defined in Table 2.on-spherical particles in medical flow problems 25 center of mass. The non-coincidence of the center of pressure and center ofmass causes the sustained oscillations (see Figure 14). Moreover, it is observedthat every particle is also slowly moving up towards the horizontal axis ofsymmetry of the domain (all particles start below the axis, see Figure 8). InFigure 15 evolution in time of the aggregated lift of five studied particles isplotted together with the evolution in time of their lift coefficients. One caneasily see that the oscillatory motion shown in Figure 14 is a direct consequenceof the lift force acting on the particles, while upward motion results from theaggregated lift being positive all the time. The behaviour of the y -velocity isalso 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. Thoselocal maxima and minima appear for angle of incidence ψ ∈ { π , π , π , π } ,so when particles are inclined to the flow direction. Particle 3, the thinnestone, is subjected to the highest deceleration, whereas particles 2 and 4 movemore smoothly. − − · − Time L i f t c o e ffi c i e n t a nd agg r e ga t e d li f t particle 1particle 2particle 3particle 4particle 5 Fig. 15
Comparison of the lift coefficient for the different particles. In bold lines we showthe aggregated lift over time. − − . . · − Time V e l o c i t y i n y d i r ec t i o n particle 1particle 2particle 3particle 4particle 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] orparticle dynamics libraries such as LAMMPS [28].All computations have been carried out on a two-socket system with IntelXeon E5-2699A v4 processors running at 2.40 Ghz.We will describe a prototypical blood-flow configuration and discuss thescaling of the implementation with respect to the number of cores. In partic-ular we will investigate the relation between computational effort used in theparticle model and in the Navier-Stokes solver. As in Section 5.1 all the param-eters have been chosen so as to resemble vessel, blood and platelet properties. The finite element model is implemented in Gascoigne 3D [4] and outlinedin Section 3.1, the discrete systems are approximated with a Newton-Krylowsolver using geometric multigrid preconditioning. Basic finite element routinesand the linear algebra workflow is partially parallelized using OpenMP, see [10]for details. Since the mesh handling and i/o are not parallelized, substantialspeedups are only reached for complex 3d problems.The particle model is based on a regular lattice mesh that covers the com-putational domain. The lattice elements of size L h × L h contain the individualparticles. Detection of particle-particle collision is limited to those particlesthat reside in the same lattice element or that belong to directly adjacentelements. This substantially helps to reduce the computational effort whichscales quadratically with the number of particles within each lattice element.Hence, we keep L h > L h must be chosen large enoughto avoid motion of particles accross multiple elements in one time step, i.e. k pd V p ≤ L h , where k pd > V p the maximum velocity of the particles. Further, the lattice mesh is basis on-spherical particles in medical flow problems 27 for parallelization since we can guarantee that no interaction between latticeelements which are separated by a complete layer can take place.
12 mm mm . mm Γ in Γ wall Γ out xy Fig. 17
Geometry of the numerical examples.
We run simulations for the 2D flow in a channel with a local narrowing of25% which should mimic a stenosed region of a blood channel. Figure 17displays 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 drivenby a Dirichlet profile on the inflow boundary Γ in given by v in ( x, y, t ) = 2 − y · v in with v in = 5 mm / s . On the wall boundary Γ wall we prescribe no-slip boundaryconditions v = 0 and on the outflow boundary Γ out we use the do-nothingoutflow condition µ∂ n v − p n = 0, see [14].The fluid viscosity is set to µ = 3 mg / mm · s. Particle and fluid densitiesare equal ρ = ρ p = 1 . / mm . Due to the fact that platelets constituteless than 1% of the blood volume [19] and the size of the domain we performsimulations with 165 000 particles.In all numerical examples the temporal step size for solving Navier-Stokesequations is k ns = 0 .
005 s, while time sub-step for particle advection is k pd =0 . P = ( L x , L y , L z , α top , α bot ) by means of thelimits indicated in (7) such that the full variety of dimensions and shape ispresent. Details on the procedure for parametrization of the particles are givenin Section 4.1. Then, 10 iterations of the interaction loop shown in Algorithm 2are performed. Hence, 10 time steps of the Navier-Stokes problem and 200 par-ticle dynamics substeps are performed. Fig. 18 shows the runtime for all 200 Particle dynamicsNavier-StokesOverallnumber of cores s p ee dup s ec o nd s Fig. 18
Left: runtime (in seconds) for the coupled Navier-Stokes particle dynamics simula-tion for an increasing number of cores. Right: parallel speedup for the complete simulationand for the Navier-Stokes finite element simulation and the particle dynamics simulationseparately. iterations. Furthermore we indicate the parallel speedup. These results showthat the allocation of computational time to the Navier-Stokes finite elementsolver and the particle dynamics system is rather balanced. While it is non-trivial to get a reasonable parallel speedup for highly efficient multigrid basedfinite element simulations (at least for simple 2d problems like this Navier-Stokes testcase), the scaling of the particle dynamics system is superior. Theseresults demonstrate that the number of particles is not the limiting factor forsuch coupled simulations.The key feature of our coupled Navier-Stokes particle dynamics schemeis the prediction of the hemodynamical coefficients by means of the previ-ously trained neural network instead of using analytical models, which arenot available, or running resolved simulations, which is not feasible for such alarge number of particles. In Fig. 19 we give details on the computational timespend in the different parts of the particle dynamics system. Besides advectionof the particles, the evaluation of the neural network is dominant, alghouthwe update the coefficients in every 10th step only. Here, a more systematicstudy of the impact of the update frequency should be performed. A furtheracceleration of the neural network evaluation is possible by using the GPUimplementation of PyTorch.The transition to more realistic 3d problems will substantially increase theeffort in both parts, finite elements and particle dynamics. For the Navier-Stokes simulation it has been demonstrated that realistic 3d blood flow situa-tions can be handled in reasonable time, see [8,9,10]. If the number of particlesis to be substantially increased, the neural network coupling for estimatinghemodynamic coefficients should be realized in a high performance packagesuch as LAMMPS [28] that allows for an efficient GPU implementation. on-spherical particles in medical flow problems 29 neural networkcollisionadvectionnumber of cores s p ee dup s ec o nd s 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 theneural network access for predicting the hemodynamical coefficients.
Suspensions of arbitrarily-shaped particles in a fluid are of great importanceboth in engineering and medical applications. However, the interaction of thenon-spherical particles with a fluid flow is a complex phenomenon, even forregularly-shaped particles in the simple fluid flows. The main difficulty lies indetermining hydrodynamic forces experienced by a particle due to their strongdependence on both particle shape and its orientation with respect to the fluidflow.In this paper, a model is successfully derived to simulate the motion ofnon-spherical particles in a non-uniform flow field, including translation androtation aspects. The model is designed to reflect platelets in a blood flow,both in terms of particle parameters and fluid configuration. The very goodagreement of these results obtained by the coupled finite element / neuralnetwork / particle dynamics simulation with state of the art documentationsin literature indicates an effectiveness of the presented approach and hencean encouraging potential toward medical applications. Furthermore, the bigimprovement over usual analytical interaction models is clearly seen as theneural network based model holds for a broad range of different shapes atany orientation. Moreover, using neural network to identify the transmissionof forces from fluid to the particles provides a possibility to adopt the modelto any desired shape of particle, making this method very promising.We have further documented details on the scaling of the approach to manyparticles, which, in 2d blood flow simplifications, matches the typical particledensity found for thrombocytes in blood flows. The computational effort is wellbalanced into the Navier-Stokes finite element part, the particle advection andthe evaluation of the neural network.
Acknowledgements
The authors acknowledge the financial support by the Deutsche Forschungs-gemeinschaft (DFG, German Research Foundation) - 314838170, GRK 2297 MathCoRe. TRacknowledges the support by the Federal Ministry of Education and Research of Germany,grant number 05M16NMA.
Conflict of interest
The authors declare that they have no conflict of interest.
References
1. Becker, R., Braack, M.: Multigrid techniques for finite elements on locally refinedmeshes. Numerical Linear Algebra with Applications , 363–379 (2000). Special Is-sue2. Becker, R., Braack, M.: A finite element pressure gradient stabilization for the Stokesequations based on local projections. Calcolo (4), 173–199 (2001)3. Becker, R., Braack, M.: A two-level stabilization scheme for the Navier-Stokes equa-tions. In: e.a. M. Feistauer (ed.) Numerical Mathematics and Advanced Applications,ENUMATH 2003, pp. 123–130. Springer (2004)4. Becker, R., Braack, M., Meidner, D., Richter, T., Vexler, B.: The finite element toolkit Gascoigne . URL
5. Braack, M., Burman, E., John, V., Lube:, G.: Stabilized finite element methods for thegeneralized oseen problem , 853–866 (2007)6. Chhabra, R., Agarwal, L., Sinha, N.: Drag on non-spherical particles: an evaluation ofavailable methods. Powder Technology (3), 288 – 295 (1999)7. Clift, R., Grace, J., Weber, M.: Bubbles, Drops, and Particles. New York; London:Academic Press (1978)8. Failer, L., Minakowski, P., Richter, T.: On the impact of fluid structure interactionin blood flow simulations,. Vietnam Journal of Mathematics (accepted, 2020). URL https://arxiv.org/abs/2003.05214
9. Failer, L., Richter, T.: A newton multigrid framework for optimal control of fluid-structure interactions. Optimization and Engineering (2020)10. Failer, L., Richter, T.: A parallel newton multigrid framework for monolithic fluid-structure interactions. Journal of Scientific Computing (2) (2020)11. Fasano, A., Sequeira, A.: Blood Coagulation. Springer International Publishing, Cham(2017)12. Filipovic, N., Kojic, M., Tsuda, A.: Modelling thrombosis using dissipative particledynamics method. Philosophical transactions. Series A, Mathematical, physical, andengineering sciences , 3265–79 (2008)13. Fogelson, A.L., Guy, R.D.: Immersed-boundary-type models of intravascular plateletaggregation. Computer Methods in Applied Mechanics and Engineering (25), 2087– 2104 (2008). Immersed Boundary Method and Its Extensions14. Heywood, J., Rannacher, R.: Finite element approximation of the nonstationary Navier-Stokes problem. iv. error analysis for second-order time discretization (3), 353–384(1990)15. Heywood, J., Rannacher, R., Turek, S.: Artificial boundaries and flux and pressureconditions for the incompressible navier-stokes equations. International Journal forNumerical Methods in Fluids (5), 325–352 (1996)16. Hlzer, A., Sommerfeld, M.: New simple correlation formula for the drag coefficient ofnon-spherical particles. Powder Technology (3), 361 – 365 (2008)17. Hlzer, A., Sommerfeld, M.: Lattice boltzmann simulations to determine drag, lift andtorque acting on non-spherical particles. Computers & Fluids (3), 572 – 589 (2009)18. Kimmritz, M., Richter, T.: Parallel multigrid method for finite element simulationsof complex flow problems on locally refined meshes. Numerical Linear Algebra withApplications (4), 615–636 (2010)on-spherical particles in medical flow problems 3119. L. Fogelson, A., Neeves, K.: Fluid mechanics of blood clot formation. Annual Reviewof Fluid Mechanics , 377–403 (2015)20. Lu, G., Third, J., M¨uller, C.: Discrete element models for non-spherical particle systems:From theoretical developments to applications. Chemical Engineering Science , 425– 465 (2015)21. Luding, S.: Collisions & Contacts between Two Particles, pp. 285–304. Springer Nether-lands, Dordrecht (1998)22. Mand, M., Rosendahl, L.: On the motion of non-spherical particles at high reynoldsnumber. Powder Technology (1), 1 – 13 (2010)23. Mody, N.A., King, M.R.: Three-dimensional simulations of a platelet-shaped spheroidnear a wall in shear flow. Physics of Fluids (11), 113302-113302-12 (2005)24. Muth, B., M¨uller, M.K., M¨uller, M., Eberhard, P., Luding, S.: Collision detection andadministration methods for many particles with different sizes. In: P. Ceary (ed.) Dis-crete Element Methods, DEM 07, 4th International Conference on Discrete ElementMethods, DEM 4 - Brisbane, Australia, pp. 1–18. Minerals Engineering Int. (2007)25. Ouchene, R., Khalij, M., Arcen, B., Tanire, A.: A new set of correlations of drag, liftand torque coefficients for non-spherical particles and large reynolds numbers. PowderTechnology , 33 – 43 (2016)26. Ouchene, R., Khalij, M., Tanire, A., Arcen, B.: Drag, lift and torque coefficients forellipsoidal particles: From low to moderate particle reynolds numbers. Computers &Fluids , 53 – 64 (2015). Small scale simulation of multiphase flows27. Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T.,Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang, E., DeVito, Z.,Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., Chintala, S.:Pytorch: An imperative style, high-performance deep learning library. In: H. Wallach,H. Larochelle, A. Beygelzimer, F. d’Alch´e Buc, E. Fox, R. Garnett (eds.) Advancesin Neural Information Processing Systems 32, pp. 8024–8035. Curran Associates, Inc.(2019)28. Plimpton, S.: Fast parallel algorithms for short-range molecular dynamics. J. Comp.Phys. , 1–19 (1995)29. Richter, T.: Fluid-structure Interactions. Models, Analysis and Finite Elements, Lecturenotes in computational science and engineering , vol. 118. Springer (2017)30. Rosendahl, L.: Using a multi-parameter particle shape description to predict the motionof non-spherical particle shapes in swirling flow. Applied Mathematical Modelling (1),11 – 25 (2000)31. Saad, Y.: Iterative Methods for Sparse Linear Systems. PWS Publishing Company(1996)32. Sanjeevi, S., Kuipers, J., Padding, J.: Drag, lift and torque correlations for non-sphericalparticles from stokes limit to high reynolds numbers. International Journal of MultiphaseFlow , 325–337 (2018)33. Shardt, O., Derksen, J.: Direct simulations of dense suspensions of non-spherical parti-cles. International Journal of Multiphase Flow , 25 – 36 (2012)34. Sweet, C.R., Chatterjee, S., Xu, Z., Bisordi, K., Rosen, E.D., Alber, M.: Modellingplatelet-blood flow interaction using the subcellular element langevin method. Journalof The Royal Society Interface (65), 1760–1771 (2011)35. Tosenberger, A., Ataullakhanov, F., Bessonov, N., Panteleev, M., Tokarev, A., Volpert,V.: Modelling of thrombus growth in flow with a dpd-pde method. Journal of TheoreticalBiology , 30 – 41 (2013)36. Viallat, A., Abkarian, M.: Dynamics of Blood Cell Suspensions in Microflows. CRCPress, Boca Raton (2019)37. Wadell, H.: The coefficient of resistance as a function of reynolds number for solids ofvarious shapes. Journal of the Franklin Institute (4), 459 – 490 (1934)38. Wiwanitkit, V.: Plateletcrit, mean platelet volume, platelet distribution width: Its ex-pected values and correlation with parallel red blood cell parameters. Clinical andApplied Thrombosis/Hemostasis (2), 175–178 (2004)39. Xu, Z., Chen, N., Kamocka, M., Rosen, E., Alber, M.: A multiscale model of thrombusdevelopment. Journal of the Royal Society, Interface / the Royal Society , 705–22(2008)2 Martyna Minakowska et al.40. Xu, Z., Chen, N., Shadden, S.C., Marsden, J.E., Kamocka, M.M., Rosen, E.D., Alber,M.: Study of blood flow impact on growth of thrombi using a multiscale model. SoftMatter , 769–779 (2009)41. Xu, Z., Lioi, J., Mu, J., Kamocka, M.M., Liu, X., Chen, D.Z., Rosen, E.D., Alber, M.: Amultiscale model of venous thrombus formation with surface-mediated control of bloodcoagulation cascade. Biophysical Journal (9), 1723 – 1732 (2010)42. Yan, S., He, Y., Tang, T., Wang, T.: Drag coefficient prediction for non-spherical parti-cles in dense gassolid two-phase flow using artificial neural network. Powder Technology , 115–124 (2019)43. Yow, H., Pitt, M., Salman, A.: Drag correlations for particles of regular shape. AdvancedPowder Technology (4), 363 – 372 (2005)44. Zastawny, M., Mallouppas, G., Zhao, F., van Wachem, B.: Derivation of drag and liftforce and torque coefficients for non-spherical particles in flows. International Journalof Multiphase Flow , 227 – 239 (2012)45. Zhang, P., Gao, C., Zhang, N., Slepian, M., Deng, Y., Bluestein, D.: Multiscale particle-based modeling of flowing platelets in blood plasma using dissipative particle dynamicsand coarse grained molecular dynamics. Cellular and Molecular Bioengineering (4)(2014)46. Zhong, W., Yu, A., Liu, X., Tong, Z., Zhang, H.: Dem/cfd-dem modelling of non-spherical particulate systems: Theoretical developments and applications. Powder Tech-nology , 108 – 152 (2016)47. Zhu, H., Zhou, Z., Yang, R., Yu, A.: Discrete particle simulation of particulate systems:Theoretical developments. Chemical Engineering Science (13), 3378 – 3396 (2007)48. Zhu, H., Zhou, Z., Yang, R., Yu, A.: Discrete particle simulation of particulate systems:A review of major applications and findings. Chemical Engineering Science63