Using a deep neural network to predict the motion of under-resolved triangular rigid bodies in an incompressible flow
UUsing a deep neural network to predict the motion of under-resolvedtriangular rigid bodies in an incompressible flow
Henry von Wahl a, ∗ and Thomas Richter a a Institute for Analysis and Numerics, Otto-von-Guericke Universität, Universitätsplatz 2, 39106 Magdeburg, Germany ∗ Corresponding author: [email protected]
Abstract
We consider non-spherical rigid body particles in an incompressible fluid in the regime where theparticles are too large to assume that they are simply transported with the fluid without back-couplingand where the particles are also too small to make fully resolved direct numerical simulations feasible. Un-fitted finite element methods with ghost-penalty stabilisation are well suited to fluid-structure-interactionproblems as posed by this setting, due to the flexible and accurate geometry handling and for allowingtopology changes in the geometry. In the computationally under resolved setting posed here, accuratecomputations of the forces by their boundary integral formulation are not viable. Furthermore, analyt-ical laws are not available due to the shape of the particles. However, accurate values of the forces areessential for realistic motion of the particles. To obtain these forces accurately, we train an artificial deepneural network using data from prototypical resolved simulations. This network is then able to predictthe force values based on information which can be obtained accurately in an under-resolved setting. Asa result, we obtain forces which are on average an order of magnitude more accurate compared to thedirect boundary-integral computation from the Navier-Stokes solution.
The aim of this work is to simulate multiple small triangular rigid bodies inside an incompressible fluidusing a moving domain approach [LO19; BFM19; WRL20], within an unfitted finite element method knownas CutFEM [Bur+14]. Solids in flows are the subject of various applications in engineering and medicine.The homogenised behaviour of these particles is well studied, but numerical simulations on the fine scale,considering resolved particles are rare, especially when non-spherical bodies are considered.If the solid particles within the fluid are very small, they can be considered as point-like within the flowfield and the interaction can be represented in a discrete element framework in terms of analytical lawsfor the fluid mechanical coefficients, as well as by modelling the effect of the bodies on the surroundingflow. Approximative laws exist for some simple shaped non-spherical particles [Fel94; HS08; MR10; Zas+12;Zho+16] and simulations with a very large number of particles can be realised in mixed continuum mechanics/ discrete element frameworks.However, if the particles are assumed to be large, only a resolved simulation remains. When consideringspherical particles, ALE approaches can be used, which allows for very good accuracy on rather coarsegrids [WT07; Wah+21] since the particle’s boundaries are exactly resolved. This approach is highly efficientand very well established in the are of elastic fluid structure interactions [HLZ81; Don82; RW10]. However,if a very large number of particles are considered, or if they are no longer radially symmetric, alternativediscretisation methods must be considered as tracking the single particles with a finite element mesh mayrequire a prohibitive computational effort. A prominent example is the immersed boundary method, goingback to Peskin [Pes02], which is able to efficiently simulate dense particle suspensions [AEW18]. Although a r X i v : . [ phy s i c s . f l u - dyn ] F e b Introduction the particles must not be exactly resolved, these methods still require an increased mesh resolution in thevicinity of the particles.Finally, we note that on a theoretical level, the motion of rigid bodies of arbitrary shapes has been studiedby Brenner and summarised by Happel and Brenner in [HB81]. However, these considerations only take intoaccount the case of creeping flows governed by the liner Stokes equations, rather than the full non-linearNavier-Stokes equations considered here.In this work we will define a hybrid finite element / (deep) neural network framework that allows toaccurately describe the interaction on under-resolved meshes. In order to write the transfer of forces betweenflow and particles with sufficient accuracy, we will use a neural network that represents the forces actingon a particle as a function of the particular flow situation instead of directly evaluating the forces from theNavier-Stokes solution. This network is trained in a presented offline phase based on prototypical resolvedsimulations. A similar approach with non-spherical but very small particles in the linear Stokes regimehas recently been described by Minakowska et al. [MRS21]. In principle this procedure can be transferredto all interface capturing schemes that allow for a robust under-resolved representation of structures in abackground mesh. Besides different CutFEM or extended finite element approaches [HH02; FD07; BE09;Bur+14] locally adjusted or parametric finite element methods [LY15; Fan13; FR14] are also an option.
Example 1 (Motivation of approach).
To illustrate that such an approach can make sense, let usconsider the following example. We compute the Schäfer-Turek 2D-1 benchmark [ST96] once using a fittedapproach together with the Babuška-Miller trick [BM84] to evaluate the drag and lift functional, and onceusing a CutFEM approach [BH14; Mas+14] together with an isoparametric mapping approach [Leh16; Leh17]to obtain the necessary geometry approximation of the discrete level set domain but with the direct evaluationof the boundary integral to realise the force values. For both the fitted and unfitted simulations we use inf-supstable Taylor-Hood elements P k / P k − , which we shall abbreviate as TH k . For the present computations wetake the order k = 3. To make the comparison as fair as possible, we consider unstructured meshes with thesame mesh size in each computation.These and all following numerical finite element computations are implemented using the finite ele-ment library NGSolve / netgen [Sch14; Sch97]. For CutFEM computations, we additionally use the add-on ngsxfem [xfem] for unfitted finite element functionality within the NGSolve framework.In Figure 1.1 we can see the velocity solution to these computations together with the computationalmeshes. Note that is is nearly impossible to distinguish the two solutions visually. In Table 1.1 we thefurther see the benchmark quantities resulting from the two computations. Here we immediately see thatthat while the values from the fitted computations are reasonably for such coarse mesh, the values resultingfrom the CutFEM computations even result in the wrong sign for the lift coefficient.Since the solution in the bulk of the domains look so indistinguishable, it is natural to ask, whetherthere are other quantities based on the solution near the obstacle, which we can compute accurately in bothsettings. To this end, we compute the average velocity in a strip of width 0 .
05 around the obstacle. Theresults from this are presented in Table 1.2. Here we see that while the fitted FEM solution resulted in valuescloser to the reference values compared to CutFEM, the difference between the accuracy of the fitted andunfitted computations are significantly smaller. Furthermore we see that we keep multiple significant figuresof accuracy in the functional value, even at values of order 10 − .We conclude from this, that while it is difficult to obtain accurate forces from the boundary integralformulation in a coarse the CutFEM computation, other features of the solution can be obtained much moreaccurately even on such coarse meshes. Therefore, if we are able to construct a mapping from flow featuresnear the obstacle to the forces acting on the obstacle, we should be able to get more accurate force values inthe under resolved CutFEM setting. (cid:78) The remainder of this text is structured as follows. In the following section 2 we describe the basic equationsgoverning the system of a rigid body in an incompressible fluid. In subsection 2.1 we then discuss the relevantphysical parameters we will consider. Section 3 covers the details of the neural network, that is the design,the generation of the draining data and the results of the training. Furthermore, we evaluate the accuracy ofthe network prediction in comparison to the direct evaluation from the Navier-Stokes solution in the setting - - Introduction Figure 1.1:
Velocity solution to the ST-2D1 problem with h max = 0 . . Top: Fitted FEM, bottom: CutFEM. Method c drag c lift ∆ p Fitted TH TH Table 1.1:
Benchmark quantities for the ST-2D1 problem computed using fitted FEM and CutFEM on meshes with h max =0 . . Reference values obtained from (accessed on 15th Feb. 2021). Method Avg. u Avg. u Fitted TH · −3 · −5 CutFEM TH · −3 · −5 Ref. 4.293 33 · −3 · −5 Table 1.2:
Average velocity in a strip of width 0.05 around the obstacle. Reference values computed using a fitted approachwith TH elements on a mesh with h max = 0 . . - - Governing equations of the generated training data. In section 4 we will then apply the resulting network in some numericalexamples and in section 5 we give some concluding remarks. Let us consider an open bounded domain Ω = F ˙ ∪ I ˙ ∪ S ⊂ R d , d ∈ { , } , divided into a d -dimensional fluidregion F , a d -dimensional solid region S and a d − I coupling the two. Thefluid in F is governed by the incompressible Navier-Stokes equations: Find a velocity and pressure ( u , p )such that ρ F (cid:0) ∂ t u + ( u · ∇ ) u (cid:1) − div σ ( u , p ) = ρ F f div( u ) = 0 (2.1)holds in F with a given fluid density ρ F , an external body force f and the Cauchy stress tensor σ ( u , p ) = µ F ( ∇ u + ∇ u T ) − Id p where µ F = ρ F ν f is the fluid’s kinematic viscosity and ν F the dynamic viscosity. The boundary conditionswhich complete the system (2.1) will be given later. We shall consider the homogeneous form of (2.1), i.e.,we take f ≡ S = ˙ S Ni =1 B i into a finite set of distinct rigid bodies (particles) B i . Forsimplicity, we assume that they each have the same homogeneous density ρ S . Let U total i = U i + ω i × r be the motion of B i , where U i is the particle’s velocity , ω i it’s angular velocity and r = x − c B the positionvector of a point in B i relative to the body’s centroid c B . These two velocities are then governed by theNewton-Euler equations m B ∂ t U i = F + F buoyant + F gravity (2.2a) I B ∂ t ω + ω × I B ω = T + T buoyant (2.2b)with the particles mass m B = ρ F vol( B ), the force and torque exerted by the fluid on the particle F = Z ∂ B i σ ( u , p ) n d S and T = Z ∂ B i r × σ ( u , p ) d S and the particles moment of inertia tensor defined by I B = ρ S Z B i ( k r k Id − r ⊗ r ) d x . The buoyancy force and torque are F buoyant = − m F g and T buoyant = − m F r bo × g where m F = ρ S vol( B i ) is the mass of the displaced fluid and r bo is the vector from the centre of mass to thecentre of buoyancy. The centre of buoyancy is defined as the centroid of the displaced fluid volume. Since theparticles are completely submersed and both the fluid and solid have a constant density, the centre of massand centre of buoyancy coincide, such that r bo = 0 and therefore T buoyant = 0. Note that equivalently, thebuoyancy effects can be included in the system by including the forces due to gravity on the right hand sideof (2.1). However, since we do not consider pressure-robust discretisations [Joh+17] here, including buoyancyin the solid and ignoring gravity on the fluid is more accurate on the discrete level. Finally, the pull due togravity is F gravity = m S g . - - Neural Network To complete the system (2.1) we need to impose boundary conditions. We shall consider Ω as a closedaquarium, and therefore take the no-slip boundary conditions u = on ∂ Ω and u = U total i on ∂ B i (2.3)which couples the fluid and solid equations. The complete system is then defined as the solution to (2.1),(2.2) and (2.3). The pressure is uniquely defined by taking p ∈ L ( F ). Remark.
We note that in two spatial dimensions, the the quadratic term ω × I B ω vanishes in (2.2b) andthe moment of inertia reduces to the scalar I B = ρ S R B k r k d x . (cid:78) We choose our fluid and solid material parameters, in order to have a setting which can be given some physicalmeaning while remaining at small to moderate Reynolds numbers.The fluid and solid parameters are chosen to approximate coarse sand in a glycerol/water mixture. Wetake a mixture of 1 part water to 4 parts glycerol at a temperature of 21°C. The resulting relevant materialparameters are summarised in Table 2.1. The fluid parameters are obtained through an online calculatortool and the density of the solid is taken as the density of quartz . The ISO standard 14688-1:2017 definescoarse sand to have a particle size of 0.63 – 2.0 mm. We therefore take the side length of our triangles to be2.0 · −3 m. Parameter µ F ρ F ν F ρ S Value 8.5679 · −2 N s m −2 · kg m −3 · −5 m s −1 · kg m −3 Table 2.1:
Fluid and solid material properties of a 1:4 water-glycerine mixture and quartz
Taking Stokes’ law for the drag acting on a sphere to estimate the terminal velocity, we get v = ρ S − ρ S µ F gr ≈ .
036 which in turn gives a Reynolds number of Re = v rν F ≈ . As we have discussed above, our aim is to train a (deep) neural network which predicts the forces acting ona small triangular rigid body, based on information retrieved from the fluid solution within the bulk of thefluid near the obstacle. In principle this should be possible, for example Raissi et. al. [RYK20] where ableto reproduce the Navier-Stokes solution from images using a physics informed neural network [RPK19] orMinakowska et. al. [MRS21] used a very simple deep neural network to predict the forces acting on bloodplatelets of different shapes in a Stokes flow.
As the work in [MRS21] is relatively similar with respect to the aim we have here, we shall take this as ourreference point for the design of our network. We note that while in [MRS21] it was the aim to use a neuralnetwork to predict the forces based on the speed and shape of very small particle in a linear Stokes flowand without back-coupling of the particles onto the fluid, we want to predict the forces based on the fluidsolution near our particle which couples back to the fluid which is governed by the non-linear Navier-Stokesequations. , accessed on 25 th Sep. 2020 , accessed on 25 th Sep. 2020 , accessed on 25 th Sep. 2020 - - Neural Network in in Input Hidden layer Hidden layer Hidden layer outoutout
Output
Figure 3.1:
Illustration of a fully connected feed-forward networks with two inputs, three hidden layers consisting of 10, 6,and 6 neurones respectively and a three outputs.
The first question we have to consider is the choice of flow features to feed into the network. To capturethe force acting on the triangular rigid body, it makes sense that the features have to be in some senselocal to the rigid body. Point evaluations of the velocity and pressure near the rigid body would be onechoice. Unfortunately, while we found that this does indeed capture the necessary information needed by aneural network, the unresolved nature of the CutFEM discretisation means that we do not have a chance ofobtaining sufficiently accurate values to feed to the network at run-time.As we have seen in Example 1, an integral of the velocity components, can be obtained accurately in acoarse CutFEM computation. We therefore choose the mean velocity, relative to the triangles velocity, in acircle around the triangle as the input. Tests on how large this circle needs to be in order to have sufficientlyaccurate values in the under resolved setting, showed that is sufficient to take the radius as twice the side-length of the rigid body. Assuming that the mesh size is not coarser than the size of the triangular rigidbody, we can then safely apply the isoparametric mapping technique [Leh16; Leh17] to accurately integrateover this domain in the CutFEM setting, without distorting the three level sets describing the sides of thetriangular body.As the forces only depend on the angle of attack between the mean (integrated) flow relative to the trianglesvelocity and the orientation of the triangle, we shall consider the input in a reference configuration. For thiswe choose the bottom side of the triangle to be parallel to the x -axis. As a result, the network learns theforces resulting from any angle of attack in this reference orientation. To obtain the physical forces, wetherefore rotate the input into the reference configuration and rotate the drag/lift predictions back into thephysical orientation (as a scalar quantity in 2d, the torque is invariant with respect to rotation). We shall consider a fully connected feed-forward network with at least three hidden layers and the ReLUactivation function, i.e., f ( x ) = max { , x } . An example of such a network can be seen in Figure 3.1. Weshall refer to networks by the number of neurones. For example, with l/m/n we denote a network with threehidden layers consisting of l , m and n neurones respectively. The number of neurones per layer and thenumber of layers we shall need for our network will be determined experimentally, by inspecting the resultsachieved during training. See subsection 3.3 below. In order to create a network, we need to generate the appropriate input and target data with which we trainthe network. - - Neural Network Neural networks are generally only accurate if they are applied in scenarios which lie in the data range withwhich the network was trained. We therefore need some information on how fast a triangular rigid body willfall under the acceleration due to gravity as described in section 2 with the material parameters described insubsection 2.1.In subsection 2.1 we calculated an estimate of the particle’s maximal velocity. To get a better sense ofthe speeds with which we need to generate the training data with, we consider an ALE discretisation of asingle triangular rigid body falling without rotation in a viscous fluid. To this end we consider the domainΩ = (0 , . × (0 , .
2) and a triangular particle of side-length 0 .
002 with the centroid (0 . , .
15) at t = 0. Forthe fluid and solid parameters we take the parameters described in Table 2.1. We then solve a simplified formof the system (2.1), (2.2) and (2.3) until t = 10. To make the ALE mapping simple, we only allow motion inthe vertical direction, i.e., we ignore the effects of horizontal drag and torque. For the construction of such amapping we refer to [Wah+21]. We consider this set up for 10 different angles of attack of the triangle andevaluate the terminal velocity.For this computation we use a mesh with h max = 1 . · − in the bulk and h max = 2 . · − on the triangleboundary with P / P dc finite elements, resulting in a FE space with dof = 2 . · and dof cond = 8 . · .The time-derivative and mesh velocity are approximated using the BDF1 scheme and the time-step used is∆ t = / .In the resulting computations, we observed that the maximal terminal velocity of the triangle was 0 . − .This gives us a good indication for the necessary training data range. In order to create a network of which we can expect sufficient accuracy, we generate the training data in asetting which is as close to the application as possible. To order to obtain accurate values for the drag, liftand torque in the case of pure translational movement, we generate the training data using a fitted ALEdiscretisation. We take the domain (0 , . with the equilateral triangular obstacle located at (0 . , .
25) inthe reference configuration. This rigid body then moved from (0 . , .
25) to (0 . , .
25) and back again overa time t end . To get a wide range of relative velocities between the triangle and the mean flow around thetriangle, we accelerate the particle at different rates, i.e., consider different values for t end . The physicallocation of the body is then given by (0 . − .
15 cos (cid:0) πt t end (cid:1) , . . for t in [0 , t end ]. To implement this motion, we use a prescribed ALE mapping, as above. In order to generatethe data of motion in different directions, we rotate the rigid body and rotate the resulting relative velocitiesand forces back into the reference configuration. As a result, we can reuse the same ALE mapping to simulatedifferent directions of relative motion of the solid body.To generate the data set, we consider t end ∈ { , . , , , , , } and rotate the triangle with angle α ∈{ iπ · | i = 0 , . . . , } . Since the triangle is equilateral, the remaining angles of attack α ∈ [ π / , π ) can beobtained by post-processing the data appropriately. For the discretisation we take TH elements on a meshwith h max = 0 .
01 and h = 1 . · − on the boundary of the rigid body. In time we discretise using a simpleBDF1 scheme with the time-step ∆ t = / . The resulting data set then contains 2 . · input/outputpairs. Remark.
It is known, that in general, the (fully) implicit BDF1 method is not a good choice to discretisethe time-derivative in the Navier-Stokes equations, since the scheme ist too diffusive, thus preventing forexample vortex shedding. However, the choice of BDF1 here is not problematic, since we are consideringflows with small Reynolds numbers. Furthermore, we have tested computations with the less diffusive BDF2scheme and no significant differences could be observed in the solution. (cid:78) dof denotes the degrees of freedom of the full finite element space. dof cond denotes the number of degrees of freedom remainingafter static condensation of unknowns internal to each element, i.e, the size of the system to be solved in every time-step. - - Neural Network Since the above ALE computations only include translational motion, this is equivalent to only consideringa parallel flow around a fixed obstacle. For the network to be universally usable, we therefore also need toinclude rotational flow data. This corresponds to rotating the triangle. However, in this situation, usingALE is not suitable, as relatively small rotations will lead to mesh-entanglement. We therefore use a highlyresolved moving domain CutFEM simulation, to generate data with rotational input. Further details of thismethod are given below in section 4To this end, we consider a triangular rigid body located at the centre of the domain Ω = (0 , . . Thisis then rotated at different speeds, clockwise and anti-clockwise, such that the total rotation at time t with respect to the initial configuration is given by sin(2 π · t/t end ). To obtain different speeds we choose t end ∈ { . , , , , , } . The mesh is constructed by taking h max = 0 . t = / together withBDF2 time-stepping. As a result, we then obtain an additional 8 . · data sets. We implement the neural network described in subsection 3.1 using
PyTorch [Pas+19]. We use the meansquared error as the loss function and take the Adam algorithm as the optimiser with a step size of 10 − .The network is trained for a total of 20000 epochs. For the network to be able to predict all three values atthe same time, we scale both the input and output data to be in the interval [ − , t end This then consists of 4 . · data points. During training we then alsoevaluate the network on the validation data set and ensure that the the training error does not decrease whilethe validation error increases.The errors of the predictions made by the networks on the training and validation data sets can be seen inTable 3.1 for the case where a single network was used to predict drag lift and torque, while Table 3.2 showsthe error for separate networks for each force. To make the results on data sets of different sizes comparable,we use the norm k err k mean := N P Ni =1 (prediction i − value i ) .To find an appropriate network size, which is large enough to capture all the information contained inthe data set, while also being small enough for fast evaluations in the final solver, we consider a number ofdifferent networks. The chosen number of layers and neurones per layer can be seen in the first column ofTable 3.1.Looking at the results in Table 3.1, we see that the results are broadly similar for all six considered networkarchitectures. Looking at the three layer networks, we observe that the prediction error resulting from the30/20/10 network are approximately 2 times larger than that of the 50/20/20 network. Looking at the fourlayer networks, we see that the errors are generally the same as for the 50/20/20 network, with some smalldeviations in both directions.In order to check whether we can gain more accuracy by considering separate networks for the threefunctionals, we train three separate 50/20/20 networks. The results thereof can be seen in Table 3.2, wherewe observe that the prediction errors are about the same as those realised by the single network with thesame architecture in Table 3.1.Both in Table 3.1 and Table 3.2 we see that the order of magnitude of the errors in both norms are verydifferent. To give these errors more meaning, we plot the predictions and the target values for a randomselection of point from the validation data set in Figure 3.2. Here we see that in general the predicted valuesmatch the target values. This that the order of magnitude of the errors in Table 3.1 and Table 3.2 are dueto the scaling of the different functional values.As a result of the above considerations, we choose the single 50/20/20 network for our application ofpredicting the drag lift and torque from the average velocity around the rigid particle. Due to the simpleinput for the network network, requiring two integrals over a relatively small area and a rotation of thesetwo values into the reference configuration, as well as the small size of the network itself, the additionalcomputation effort introduced to the solver by predicting the forces, rather than evaluating them directly, is - - Neural Network Training data Validation dataArchitecture Unknowns Functional k · k mean k · k ∞ k · k mean k · k ∞ F · −3 · −2 · −3 · −2 Force F · −3 · −2 · −3 · −2 Torque 1.5 · −6 · −5 · −6 · −5 F · −3 · −2 · −3 · −2 Force F · −3 · −2 · −3 · −2 Torque 1.5 · −6 · −5 · −6 · −5 F · −3 · −2 · −3 · −2 Force F · −3 · −2 · −3 · −2 Torque 1.5 · −6 · −5 · −6 · −5 F · −3 · −2 · −3 · −2 Force F · −3 · −2 · −3 · −2 Torque 1.5 · −6 · −5 · −6 · −5 F · −3 · −2 · −3 · −2 Force F · −3 · −2 · −3 · −2 Torque 1.5 · −6 · −5 · −6 · −5 Table 3.1:
Prediction errors in a wighted ‘ norm and the maximum norm on the training and validation data sets after20000 epochs of training for a number of different network sizes. Each network predicts all three functionalvalues simultaneously. Training data Validation dataFunctional k · k mean k · k ∞ k · k mean k · k ∞ Force F · −3 · −2 · −3 · −2 Force F · −3 · −2 · −3 · −2 Torque 1.5 · −6 · −5 · −6 · −6 Table 3.2:
Prediction errors in a wighted ‘ norm and the maximum norm on the training and validation data sets after20000. Separate networks of the architecture 50/20/20 with 1611 unknowns are used top predict the threefunctional values. − · − Force F − · − Force F − · − TorqueTarget data Prediction
Figure 3.2:
Target data and network prediction for 300 random points in the training data set. Network architecture:50/20/20, 1653 unknowns. - - Numerical examples Drag Lift TorqueForce comp. Discr. k · k mean k · k ∞ k · k mean k · k ∞ k · k mean k · k ∞ Evaluation TH · −2 · −1 · −3 · −2 · −4 · −3 Prediction TH · −3 · −2 · −4 · −2 · −7 · −6 Evaluation TH · −2 · −1 · −3 · −1 · −4 · −2 Prediction TH · −3 · −2 · −4 · −2 · −7 · −6 Table 3.3:
Absolute errors in the forces resulting from the direct evaluation of the boundary integrals and the predictionmade by the neural network in a CutFEM computation. negligible compared to the effort required to solve the no-linear system resulting from the FEM discretisationof the Navier-Stokes equations.
Remark (Training times).
The training the above networks was performed on a Tesla V100 PCIe 16GBgraphics card with
PyTorch using
CUDA version 10.1. Due to the small size and simple structure of thenetworks, the training times for 2.0 · epochs ranged between 277 seconds for the 30/20/10 network and936 seconds for the 100/50/20/20 network. (cid:78) In the previous section we have compared the network predictions against the “true” values in the sense of theALE training data. However, the aim for the network is to predict the force values in a CutFEM simulationand to be more accurate than the boundary integral evaluation.To test the accuracy of the network in the CutFEM setting, we run a moving domain CutFEM simulationof the training data generation set up for t end = 3. This is in order to have an fitted reference solution for theHere we took a background mesh with h max = 10 − in the area where we compute the average of the velocityaround the triangle and h max = 0 .
04 in the remaining part of the domain. The mesh in the averaging area istherefore a factor of 2 smaller than the size of the rigid body. On this mesh we work with unfitted TH and TH elements. We then compute the errors of the force prediction and evaluation, by comparing the valuesagainst the direct evaluation of a highly resolved ALE computation of the identical set-up. Here we took h max = 0 .
01 in the bulk and h max = 10 − on the triangle. This corresponds to 20 elements on each triangleside. On this fitted mesh we used TH elements. In both cases, the time-step ∆ t = / was chosen.The resulting prediction errors of the forces in the CutFEM simulation can be seen Table 3.3. Here we seethat the mean and maximal predictions is at least one order of magnitude smaller than the direct evaluation.In Figure 3.3 we plot the resulting forces for single run of the above computations ( α = 0). This shows, thatwhile there is a visible difference between the prediction and the “ground truth”, the predictions mirror thereal forces much better than the direct computation via the boundary integral evaluation. Here we also note,that there was no significant improvement in the direct evaluation when using higher order elements. Weconclude that while the predictions are not perfect, we have improved on the direct computation of the forces.However, we also note that we cannot expect any asymptotic mesh convergence here, as the prediction errorwill begin to dominate once the interface is sufficiently resolved in every time-step. We consider a number of numerical examples which use the neural network trained in the previous sectionin under resolved, moving domain, unfitted finite element method. The details and analysis of this methodas applied to the time-dependent Stokes equations on moving domain is given in [WRL20]. The idea of thismethod is use ghost-penalty stabilisation to implement a discrete, implicit extension of the solution into theexterior of the fluid domain, in order to apply a simple BDF formula to the time derivative in the case ofa moving interface. The only modification here is the addition of the convective term. See also [Wah+21], - - Numerical examples − . − . . . Force F − − − · − Force F − − − · − TorqueALE - Evaluation CutFEM - Evaluation CutFEM - Network prediction
Figure 3.3:
Prediction and evaluation in a CutFEM simulation with TH elements compared against the values evaluatedin a fitted ALE computation. where this method has also been applied to a fluid-structure interaction problem with contact and in whichthe fluid is described by the incompressible Navier-Stokes equations.To solve the coupled fluid-solid system, we use a partitioned approach. That is we solve the fluid andsolid equations separately and iterate between the two until the full system is solved implicitly. For the solidvelocity we use a relaxation scheme together with Aitken’s ∆ method to accelerate the relaxation scheme.A relaxation is necessary for stability of the partitioned solver of the coupled system and different relaxationstrategies, including Aitken’s ∆ method are discussed in the literature [KW08]. Remark 2 (Geometry description in CutFEM).
In CutFEM, numerical integration on elements whichare cut by the level-set function describing the geometry, is based on a piecewise linear interpolation of the(smooth) level set function. This is necessary to robustly construct quadrature rules on cut elements, see[Bur+14, Section 5.1]. To represent the three sides of the triangular rigid bodies, we therefore use threeseparate level set functions φ , φ and φ . ngsxfem [xfem] then provides the functionality to integrate withrespect to each of these level sets, for example the domain where all level sets are negative. Since the levelsets describing the sides are linear. the numerical integration on the geometry posed here is exact and we donot have to construct a single level set function φ = max( φ , φ , φ ) which would then be interpolated intothe P space on the mesh, which in turn would remove the sharp corners of the particles. (cid:78) In subsection 3.4 we compared the resulting forces from the network against the forces computed with anALE and CutFEM simulation in the setting of a prescribed velocity of the triangular rigid body. To establishthe accuracy in the coupled setting we shall consider the scenario used in subsubsection 3.2.1 for the a prioricomputations to generate the training data, where horizontal and rotational motion was neglected, in orderto compare the motion resulting from the forces predicted by the network against an easily computablereference solution.To this end we take the domain Ω = (0 , . , assert no-slip boundary conditions at the left and rightwall and the do-nothing condition at the top and bottom. The rigid body is an equilateral triangle withside-length 2 · − . The fluid and solid material parameters are as in Table 2.1. At t = 0 the centre of mass ofthe rigid body is at (0 . , .
08) and we rotate the body by an angle α with respect to reference configuration,in which the bottom of the triangle is parallel to the x -axis. We shall consider α = 0 , π/ , π/ TH elements on a mesh with h max = 5 . · − , h = 1 . · − ina horizontal strip of height 8 . · − around the rigid body and h = 4 . · − on the interface of the rigidbody.For both CutFEM computations we consider a mesh with h max = 1 . · − on which we consider TH ele-ments. For the unfitted discretisation used, we take the Nitsche parameter to enforce the Dirichlet boundary - - Numerical examples conditions on the level set interfaces as γ N = 100 while the stability and extension ghost-penalty parametersare γ gp,stab = 0 . γ gp,ext = 0 .
001 respectively. For details on these parameters we refer to [WRL20].We take the time-step ∆ t = / and compute until t = 1 .
0. For the partitioned fluid-solid solver, wetake the initial relaxation parameter as ω = 0 .
5, allow a maximum of 10 sub-iterations per time step and atolerance of 1% in the relative update.The forces acting from the fluid onto the rigid body are evaluated using the single 50/20/20 network fordrag, lift and torque. As discussed above we use an isoparametric approach to accurately integrate the fluidvelocity relative to the solid velocity in the circle of radius 2 × body-side-length.In Figure 4.1 we can see the vertical force, the vertical speed and the vertical position of the centre of mass ofthe rigid body over time for each of the three triangle orientations. Here we see that for all three orientations,the drag obtained by the network prediction is significantly more accurate than the direct evaluation of theboundary integral. Looking at the body’s vertical speed and position, we see that there is a visible differencebetween the ALE reference solution and the movement resulting from the predicted forces. The uniformlyfaster speed in the CutFEM prediction computations can be attributed to the fact, that at small speeds ofthe particle, the network underestimated the drag, as can be seen in the left column of Figure 4.1. Howeverwe also seen that the predictions are significantly better than the result from the direct evaluation of theforces. We now consider the full fluid-rigid body system including translational movement in all directions as wellas rotational motion. We begin with a single particle. The set up for this is as follows. We again take thedomain Ω = (0 , . and we assert no-slip boundary conditions at the left and right wall, while we apply ado-nothing condition at the top and bottom. This corresponds to the domain being a section of a channel.The rigid body is an equilateral triangle with side-length 2 · − such that the top side is parallel to the x -axis. The initial centre of mass of the rigid body is taken as (0 . , .
08) and the initial condition for thefluid and solid are ( u , p ) = ( , U, ω ) = (0 ,
0) respectively. The fluid and solid material parametersare taken as in Table 2.1. The resulting system is then driven by the acceleration due to gravity acting onthe on the rigid body.We take a mesh of the background domain with h max = 1 . · − (half the rigid body’s side-length) andsecond mesh with h max = 7 . · − to investigate the mesh dependence. On each mesh we take TH elements.The remaining discretisation parameters are chosen as in the unfitted simulation (with forces from the neuralnetwork) as in subsection 4.1.In Figure 4.2 we see the resulting velocity components which make up the movement of the rigid body forthe two meshes considered. The fluid solution at time t = 1 . t = 0 . t = 1 . · − , we find that the resulting maximal velocity at the corners of the triangular rigid bodies isof order 10 − . So in general, the velocity of the triangle resulting from rotation is smaller than the verticalvelocity component, which we consider to be realistic.Finally, we note that a comparison to a fitted ALE simulation is out of scope here, since this would haveto include re-meshing procedures and projections of the solution onto the resulting new meshes, in order toavoid mesh-angling resulting from the rotation of the particles. - - Numerical examples . · − α = Vertical force F . − − − − · − Vertical velocity . · − Vertical position . . . · − α = π / . − − − − · − . · − . · − α = π / . − − − − · − . · − ALE - Evaluation CutFEM - Evaluation CutFEM - Prediction
Figure 4.1:
Vertical force, velocity and position of a single falling triangular rigid body restricted to vertical motion.Comparison between an ALE reference computation, the force computation in an under resolved CutFEMcomputation and the network prediction in an under resolved CutFEM computation. - - Numerical examples . . − · − Velocity x -component . . − − · − Velocity y -component . . − Angular velocity h max = 10 − h max = 7 . · − Figure 4.2:
Translational and rotational velocity components of a single triangular rigid body in free fall on an under-resolvedmesh with the forces from the predictions by a deep neural network.
Figure 4.3:
Velocity field (left) and pressure with mesh (right) at t = 1 . resulting from a single triangular rigid body in freefall on a mesh with h max = 10 − with the forces governing the solid motion obtained from a neural network. - - Conclusions and Outlook We consider the same basic set-up as in subsection 4.2. However, we take 5 rigid triangular rigid bodies. Wedenote these by tri i , 1 = 0 , . . . ,
4. Their initial position and orientation are given by their centre of mass as theangle of rotation with respect to the reference configuration, in which the bottom side of the triangle is paral-lel to the x -axis. We denote the initial position and rotation by ( c x , c y , α ). The initial states of tri , . . . , tri are then (0 . , . , . , . , π/ . , . , π ), (0 . , . , π/
7) and (0 . , . , π/
13) respect-ively. Looking at Figure 4.5, this enumeration corresponds to starting with the top-left triangle and countingcounter-clockwise. The discretisation parameters of the moving domain CutFEM method are all as in sub-section 4.2.The results from these computations are shown in Figure 4.4 and Figure 4.5. The translational androtational velocity components are shown in Figure 4.4 while the fluid solution at t = 1 . tri has the largest (vertical) velocity. This also makes sense, as this is directly in the wakeof tri . Furthermore, we note that the angular velocity of tri and tri , which are above the other particles,is smaller than that of the other particles. We attribute this to the fact, that these two particles are in theparallel flow wake of the other three particles, resulting i a reduced torque acting on these two particles. Conclusions.
We have presented a framework using a hybrid finite element / neural network approach toobtain improved drag lift and torque values acting on a triangular rigid body. For this, we trained a smalldeep neural network using training data obtained from prototypical motion configurations in resolved finiteelement simulations. We considered a range of networks by varying the number of layers as well as the numberof neurones per layer. Here we saw that the mapping from the average velocity around the particle to theforces acing on this particle was accurately captured by relatively small networks. Applying one of thesenetworks in an under resolved moving domain CutFEM simulation, we showed that this approach results insignificantly more accurate force values resulting realistic motion, compared to the case where the forces arecomputed directly from the solution on an under resolved mesh.
Outlook.
A number of interesting questions remain for future research. For example, it would be veryinteresting compare the particle motion realised in subsection 4.3 to fitted ALE simulations using re-meshingtechniques.Furthermore, it remains to generalise this approach to other particle shapes and sizes. Here it would alsobe interesting to see, if a different choice of input data could generalise the approach to use a single neuralnetwork to predict accurate forces with different material parameters.
Acknowledgements
The authors acknowledge support by the Deutsche Forschungsgemeinschaft (DFG, German Research Found-ation) - 314838170, GRK 2297 MathCoRe.
References [AEW18]
M. Hazmil A. Azis , F. Evrard and
B. van Wachem . ‘An immersed boundary method for flows with dense particlesuspensions’. In:
Acta Mech. doi : . - -eferences . · − P a r t i c l e t r i Velocity x -component . − − · − Velocity y -component . − Angular velocity . − − . . · − P a r t i c l e t r i . − − · − . − . − − . · − P a r t i c l e t r i . − − − · − . − . − · − P a r t i c l e t r i . − − − · − . − . − − · − P a r t i c l e t r i . − − − · − . − h max = 10 − h max = 7 . · − Figure 4.4:
Translational and rotational velocity of five triangular rigid bodies in free fall in an under-resolved CutFEMsimulation where the forces governing the solid motion are from the evaluation of a deep neural network. - -eferences Figure 4.5:
Velocity field (left) and pressure with mesh (right) at t = 1 . resulting from five triangular rigid bodies in freefall on a mesh with h max = 10 − with the forces governing the solid motion obtained from a neural network. [BM84] I. Babuška and
A. Miller . ‘The post-processing approach in the finite element method-part 1: Calculation ofdisplacements, stresses and other higher derivatives of the displacements’. In:
Internat. J. Numer. Methods Engrg. doi : .[BE09] P. Bastian and
C. Engwer . ‘An unfitted finite element method using discontinuous Galerkin’. In:
Internat. J.Numer. Methods Engrg. doi : .[Bur+14] E. Burman , S. Claus , P. Hansbo , M. G. Larson and
A. Massing . ‘CutFEM: Discretizing geometry and partialdifferential equations’. In:
Internat. J. Numer. Methods Engrg. doi : .[BFM19] E. Burman , S. Frei and
A. Massing . Eulerian time-stepping schemes for the non-stationary Stokes equations ontime-dependent domains . 7th Oct. 2019. arXiv: .[BH14]
E. Burman and
P. Hansbo . ‘Fictitious domain methods using cut elements: III. A stabilized Nitsche method forStokes’ problem’. In:
ESAIM Math. Model. Numer. Anal. doi : .[Don82] J. Donea . ‘An arbitrary Lagrangian-Eulerian finite element method for transient dynamic fluid-structure interac-tions’. In:
Comput. Methods Appl. Mech. Engrg.
33 (1982), pp. 689–723. doi : .[Fan13] X. Fang . ‘An Isoparametric Finite Element Method for Elliptic Interface Problems with Nonhomogeneous JumpConditions’. In:
WSEAS Trans. Math.
12 (2013). issn : 2224-2880. url : (visited on 22/02/2021).[Fel94] R. Di Felice . ‘The voidage function for fluid-particle interaction systems’. In:
Int. J. Multiph. Flow
20 (1994),pp. 153–159. doi : .[FD07] K. Fidkowski and
D. Darmofal . ‘A triangular cut-cell adaptive method for high-order discretizations of thecompressible Navier–Stokes equations’. In:
J. Comput. Phys. doi : .[FR14] S. Frei and
T. Richter . ‘A locally modified parametric finite element method for interface problems’. In:
SIAMJ. Num. Anal.
52 (2014), pp. 2315–2334. doi : .[HH02] A. Hansbo and
P. Hansbo . ‘An unfitted finite element method, based on Nitsche’s method, for elliptic interfaceproblems’. In:
Comput. Methods Appl. Mech. Engrg. doi : .[HB81] J. Happel and
H. Brenner . Low Reynolds number hydrodynamics . Springer Netherlands, 1981. doi : .[HS08] A. Hölzer and
M. Sommerfeld . ‘New simple correlation formula for the drag coefficient of non-spherical particles’.In:
Powder Technol. doi : .[HLZ81] T. J. R. Hughes , W. K. Liu and
T. K. Zimmermann . ‘Lagrangian-Eulerian finite element formulations for in-compressible viscous flows’. In:
Comput. Methods Appl. Mech. Engrg.
29 (1981), pp. 329–349. doi : .[Joh+17] V. John , A. Linke , C. Merdon , M. Neilan and
L. G. Rebholz . ‘On the divergence constraint in mixed finiteelement methods for incompressible flows’. In:
SIAM Rev. doi : . - -eferences [KW08] U. Küttler and
W. A. Wall . ‘Fixed-point fluid-structure interaction solvers with dynamic relaxation’. In:
Comput.Mech. doi : .[LY15] U. Langer and
H. Yang . Numerical simulation of parabolic moving and growing interface problems using smallmesh deformation . 31st July 2015. arXiv: .[Leh16]
C. Lehrenfeld . ‘High order unfitted finite element methods on level set domains using isoparametric mappings’.In:
Comput. Methods Appl. Mech. Engrg.
300 (Mar. 2016), pp. 716–733. doi : .[Leh17] C. Lehrenfeld . ‘A higher order isoparametric fictitious domain method for level set domains’. In:
Lecture Notes inComputational Science and Engineering . Ed. by
S. Bordas , E. Burman , M. Larson and
M. Olshanskii . Vol. 121.Springer, Cham, 2017, pp. 65–92. doi : .[LO19] C. Lehrenfeld and
M. A. Olshanskii . ‘An Eulerian finite element method for PDEs in time-dependent domains’.In:
ESAIM Math. Model. Numer. Anal. doi : .[MR10] M. Mandø and
L. Rosendahl . ‘On the motion of non-spherical particles at high Reynolds number’. In:
PowderTechnol. doi : .[Mas+14] A. Massing , M. G. Larson , A. Logg and
M. E. Rognes . ‘A stabilized Nitsche fictitious domain method for theStokes problem’. In:
J. Sci. Comput. doi : .[MRS21] M. Minakowska , T. Richter and
S. Sager . ‘A finite element / neural network framework for modeling suspensionsof non-spherical particles.’ In:
Vietnam J. Math. (Feb. 2021). doi : .[xfem] ngsxfem : An add-on to NGSolve for unfitted finite element discretizations . Version 1.3.2101. 2021. url : https ://github.com/ngsxfem/ngsxfem .[Pas+19] A. Paszke et al. ‘PyTorch: An Imperative Style, High-Performance Deep Learning Library’. In:
Advances in NeuralInformation Processing Systems 32 . Ed. by
H. Wallach et al. Vol. 2. Curran Associates, Inc., 2019, pp. 8026–8037. url : http://papers.nips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf (visited on 22/02/2021).[Pes02] C. S. Peskin . ‘The immersed boundary method’. In:
Acta Numer.
11 (Jan. 2002), pp. 479–517. doi : .[RPK19] M. Raissi , P. Perdikaris and
G. E. Karniadakis . ‘Physics-informed neural networks: A deep learning frameworkfor solving forward and inverse problems involving nonlinear partial differential equations’. In:
J. Comput. Phys.
378 (Feb. 2019), pp. 686–707. doi : .[RYK20] M. Raissi , A. Yazdani and
G. E. Karniadakis . ‘Hidden fluid mechanics: Learning velocity and pressure fields fromflow visualizations’. In:
Science doi : .[RW10] T. Richter and
T. Wick . ‘Finite Elements for Fluid-Structure Interaction in ALE and Fully Eulerian coordinates’.In:
Comput. Methods Appl. Mech. Engrg. doi : .[ST96] M. Schäfer and
S. Turek . ‘Benchmark computations of laminar flow around a cylinder’. In:
Notes on NumericalFluid Mechanics (NNFM) . Vieweg+Teubner Verlag, 1996, pp. 547–566. doi : .[Sch97] J. Schöberl . ‘NETGEN an advancing front 2D/3D-mesh generator based on abstract rules’. In:
Comput. Vis. Sci. doi : .[Sch14] J. Schöberl . C++11 implementation of finite elements in NGSolve . ASC Report No. 30/2014. Tech. rep. Institutefor Analysis and Scientific Computing, TU Wien, 26th Sept. 2014. url : (visited on 22/02/2021).[Wah+21] H. von Wahl , T. Richter , S. Frei and
T. Hagemeier . ‘Falling balls in a viscous fluid with contact: Comparingnumerical simulations with experimental data’. In:
Phys. Fluids doi : .[WRL20] H. von Wahl , T. Richter and
C. Lehrenfeld . An unfitted Eulerian finite element method for the time-dependentStokes problem on moving domains . 6th Feb. 2020. arXiv: .[WT07]
D. Wan and
S. Turek . ‘Fictitious boundary and moving mesh methods for the numerical simulation of rigidparticulate flows’. In:
J. Comput. Phys. doi : .[Zas+12] M. Zastawny , G. Mallouppas , F. Zhao and
B. van Wachem . ‘Derivation of drag and lift force and torquecoefficients for non-spherical particles in flows’. In:
Int. J. Multiph. Flow
39 (Mar. 2012), pp. 227–239. doi : .[Zho+16] W. Zhong , A. Yu , X. Liu , Z. Tong and
H. Zhang . ‘DEM/CFD-DEM Modelling of non-spherical particulatesystems: Theoretical Developments and Applications’. In:
Powder Technol.
302 (Nov. 2016), pp. 108–152. doi : . -18