# Lehrstuhl für Informatik 10 (Systemsimulation) · 2.1 Rigid Body Dynamics ... 2D triangle ......

### Transcript of Lehrstuhl für Informatik 10 (Systemsimulation) · 2.1 Rigid Body Dynamics ... 2D triangle ......

FRIEDRICH-ALEXANDER-UNIVERSITÄT ERLANGEN-NÜRNBERGTECHNISCHE FAKULTÄT • DEPARTMENT INFORMATIK

Lehrstuhl für Informatik 10 (Systemsimulation)

Complex Crystalline Geometries - Simulating Particle Filters

Vasanth Baskar

Master Thesis

Complex Crystalline Geometries - Simulating Particle Filters

Vasanth BaskarMaster Thesis

Aufgabensteller: Prof. Dr. U. RüdeBetreuer: M. Sc. S. Eibl, M. Sc. C. RettingerBearbeitungszeitraum: 16.05.2017 – 16.11.2017

Erklärung:

Ich versichere, dass ich die Arbeit ohne fremde Hilfe und ohne Benutzung anderer als der angegebenenQuellen angefertigt habe und dass die Arbeit in gleicher oder ähnlicher Form noch keiner anderen Prü-fungsbehörde vorgelegen hat und von dieser als Teil einer Prüfungsleistung angenommen wurde. AlleAusführungen, die wörtlich oder sinngemäß übernommen wurden, sind als solche gekennzeichnet.

Der Universität Erlangen-Nürnberg, vertreten durch den Lehrstuhl für Systemsimulation (Informatik 10),wird für Zwecke der Forschung und Lehre ein einfaches, kostenloses, zeitlich und örtlich unbeschränktesNutzungsrecht an den Arbeitsergebnissen der Master Thesis einschließlich etwaiger Schutzrechte undUrheberrechte eingeräumt.

Erlangen, den 16. November 2017 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

AbstractSimulation of particulate flow through fibrous media is of great importance in many industries. A

fluid-solid coupling to simulate the flow of particles through the porous media has been developed tostudy the importance of the fiber characteristics and its effects on the particle deposition and its efficiency.The scope of this thesis is to implement a new convex rigid body described using triangles, representing thefilter crystals. It is distributed randomly to form a structure that resembles a porous media through whichparticles are driven. A critical step in designing this rigid body also involves implementing and validatingthe collision detection routine between crystals and spheres. An in-house developed framework known asPhysics Engine (pe) is used for the simulations involving rigid bodies. waLBerla ( Widely applicable LatticeBoltzmann from Erlangen ) is used for the fluid simulations. The coupling between pe and waLBerlaprovides with an opportunity to incorporate both the rigid bodies and fluid flow resulting in the particle-based direct numerical simulation. Influence of Solid Volume Fraction (SVF) on the filtration efficiencywill be determined.

4

Contents

Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

List of algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81.2 Overview and Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.1 Rigid Body Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.1.1 Physics Engine . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.2 Lattice Boltzmann Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.2.1 waLBerla . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162.3 Coupling - Momentum Exchange Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2.3.1 pe - waLBerla Coupling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

3 Domain setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 193.1 Crystal Body . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

3.1.1 Convex Polyhedra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 193.1.2 Creation and Import . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203.1.3 Collision Detection between Crystal Body and Sphere . . . . . . . . . . . . . . . . . 233.1.4 Mapping in waLBerla . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

4 Sizing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 284.1 Crystal Body Sizing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 284.2 Particle Sizing and Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

5 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 315.1 Results with pe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 315.2 Results with pe-waLBerla Coupling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 436.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 436.2 Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

7 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

5

List of Figures

2.1 Relative velocities : Contact Point . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.2 Collision Parameters Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142.3 D3Q19 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152.4 Representation of Objects in pe and waLBerla . . . . . . . . . . . . . . . . . . . . . . . . . . 182.5 Discretization error of Mapping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183.1 Convex Polygon . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 193.2 Edges and Face Normal of a Triangle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203.3 Texture of ACM Diesel Particulate Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203.4 Filter - Single Crystal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 213.5 Single Crystal - 3D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 213.6 Base Crystal Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 213.7 Class Dependency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223.8 Prototype - Filter Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 243.9 Voronoi Region - 2D triangle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 243.10 Mapping of Crystal to waLBerla . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 264.1 Domain size in waLBerla . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 294.2 Domain with Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 305.1 Penetration Test 1(a) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 315.2 Penetration Test 1(b) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 325.3 Penetration Test 1(c) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 325.4 Penetration Test 2(a) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 325.5 Penetration Test 2(b) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 335.6 Penetration Test 2(c) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 335.7 Collision test 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 345.8 Collision Test 2(a) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 345.9 Collision Test 2(b) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 355.10 Collision Test 2(c) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 355.11 Flag Field . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 355.12 Velocity Profile . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 365.13 Filtration Efficiency vs SVF - Plot 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 375.14 Filter with trapped particles - Setup 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 375.15 Filter with trapped particles - Setup 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 375.16 Filtration Efficiency vs SVF - Plot 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 385.17 Different Filter Sizes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 385.18 Solid Volume Fraction Change with different ∆x . . . . . . . . . . . . . . . . . . . . . . . . . 395.19 Filtraton Efficiency based on different discretization size of particles and domain . . . . . 395.20 Domain with One layer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 405.21 Domain with Two layers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 405.22 Domain with Three layers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 405.23 Filtration Efficiency for Different Layers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 415.24 Particles trapped in one layer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 415.25 Particles trapped in two layers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 425.26 Particles trapped in three layers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

6

List of Algorithms1 pe Rigid body simulation time step . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 LBM Sweep . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173 LBM Algorithm with Momentum Exchange Method . . . . . . . . . . . . . . . . . . . . . . . 184 Finding the Closest Point to Triangle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 255 Collide Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 256 Mapping Of Crystal Body to waLBerla . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

7

1 Introduction

1.1 MotivationFluid flow in porous media has a variety of applications namely air conditioning systems, purification

processes in petroleum products and also limiting the emission of particulate matter from diesel engines.A Diesel engine is an internal combustion engine widely used in the on-road transport because of its highefficiency and reliability. Fresh air is compressed and at the end of the compression cycle fuel is injected,converting chemically bound energy into useful mechanical work. During the combustion, in additionto combustion byproducts such as carbon dioxide (CO2) and water (H2O) undesirable pollutants suchas carbon monoxide (CO), hydrocarbons (HC), nitrous oxides (NOx ), and particulate matters (PM) areproduced (Fiebig et al., 2014).

The work of the Diesel Particulate Filter (DPF) is to trap the particulate matters before it entersthe atmosphere. These particulate matters comprise mainly of ash and soot (Chen et al., 2016). Thisaccumulation of ash and soot leads to the increase in back pressure and eventually drop in engineperformance (Bensaid et al., 2010). Recent studies of the filter performance have led to the developmentof fibrous filtration techniques, where irregularly shaped fibers are arranged randomly to improve theperformance as well as limit the pressure drop across the filter (Zhu et al., 2000). With many layers ofrandomly arranged fibers, the filters can approach the efficiency of 100% (Karadimos and Ocone, 2003).

Various researchers have studied the fluid-solid flow numerical simulations involving specific arrange-ments of the fibrous filter media such as Karadimos and Ocone (2003), Przekop et al. (2004), and Dunnettand Clement (2006). Novel approaches have been undertaken to understand the performance of thefilters, which can be characterized by the filter efficiency based on the permeability or porosity and alsothe fiber diameters. Müller et al. (2014) discusses the direct numerical simulation carried out for fluid fieldcalculations using CFD. Simulation of pressure drop across the filter using the ANSYS-Fluent CFD codealong with the calculation of the collection efficiency is discussed in Hosseini and Tafreshi (2012). Effect ofthe change in fiber parameters on the particle capture efficiency using the CFD-DEM model is discussedin Qian et al. (2013).

With the traditional CFD methods, it is difficult to simulate the fluid-particle interactions. This isbecause for the moving particles the fluid domain has to be continuously re-meshed increasing thedifficulty. Lattice Boltzmann Method (LBM) has been used to simulate the particulate flows because of itsflexibility in the treatment of the moving boundary conditions. A coupling between the LBM and DiscreteElement Method (DEM) to simulate the flow in the porous media has already been discussed by Qiu(2015).

The chair of System Simulation from the department of Computer Science at the University of Erlangen-Nuremberg has been developing a rigid multi-body physics engine framework called as Physics Engine(pe) since 2006. This pe is primarily used for simulation in the field of granular media where an arbitrarynumber of particles interact with each other. It provides a realistic platform for the simulation of billions ofrigid bodies in a virtual reality environment with this massively parallel framework (Iglberger, 2010). Thereare two different collision operators available in pe namely Discrete Element Method (DEM) implementedby Heene (2011) and time stepping Fast Frictional Dynamics (FFD) solver called as Hard Contact Semi-Implicit Time Stepping (HCSITS) solver implemented by Preclik (2014).

Along with pe, another in-house developed framework called waLBerla (widely applicable LatticeBoltzmann from Erlangen) is used for the coupling of the rigid bodies with the fluid. The coupling ofrigid body framework along with waLBerla is used to simulate the movement of particles under thehydrodynamic forces. This coupling using the momentum exchange methods have been studied beforefor the simulation of particulate flows (Rettinger and Rüde, 2017).

The purpose of this thesis is to create a new rigid body type in pe called the crystal body, which isthe fundamental part of the crystal structure of the diesel particulate filter. This crystal body is a convexstructure described using the vertices of the triangles which comprises the faces of this convex structure.These triangles are used in the development of the collision detection and response algorithm for thisrigid body. This developed collision operator is only for interaction between the spherical particles andthis convex crystal body. This limitation has to be understood. Then this convex structure is mapped inthe fluid-solid coupling as a fixed object with moving particles dispersed in the flow. With heterogeneousconfigurations and alignments of the crystal structure, the efficiency of the filter is studied.

8

1.2 Overview and StructureSection 2 discusses the theoretical background and the introduction to the frameworks used in this

work. The coupling, momentum exchange method, will be broached in this section . Section 3 and 4addresses how the rigid body is generated and imported in the pe. It also examines the algorithm used toassign the solid cells in the lattice structure of the fluid simulation. Section 5 validates and explains theresults. Finally, section 6 reviews the work done in this thesis and gives a possibilities of moving forward.

9

2 Background

2.1 Rigid Body DynamicsParticle and rigid body systems are the most commonly used method to simulate physically based

simulations. Particle systems consider just the point mass whereas rigid body systems consider the massdistribution and shape of the objects. Rigid body simulations have been used in diverse fields from vehicledynamics to physics engines in computer games. These simulations have been used to observe the realisticbehavior of objects under the influence of mechanical forces. The motion of rigid bodies can be classifiedas unconstrained or constrained. Unconstrained motion can be described as the movement under theinfluence of Newton’s law of motion without the consideration of the contact with other rigid bodies.Constrained motions on the other hand are characterized as such when the rigid bodies are not allowed tohave any inter-penetration between each other. The movement of rigid bodies dealing with constrained orunconstrained motion has been discussed in Baraff (2001). Some basic concepts of rigid body dynamicsare explained well in Kavan (2003).

All rigid bodies occupy a volume of space and a particular shape. It cannot undergo deformation,meaning that it can only have translational and rotational movements. Let us denote the body’s volumeV ⊆R3 and density function ρ : V 7→R, then we can define the body’s mass as

m =∫

Vρ(r )dV (2.1)

and center of mass ∫V rρ(r )dV

m(2.2)

and the center of mass relative to the body space vector can be written as

rc = (rc,x ,rc,y ,rc,z ) (2.3)

where subscript "c" refers to the vectors relative to the center of mass. At the start of each simulation step,the dynamic state of all the rigid bodies is known such as their linear and angular positions, velocities, andaccelerations. The position of the rigid body is characterized by translation vector and a rotation matrix.

The position of the rigid body can be described relative to its center of mass. At time t we define theposition of the vector with respect to the world space rw (t ) as

rw (t ) = rc (t )+R(t )rb (2.4)

where rc (t ) is the center of mass concerning the world space and rb(t ) is the vector from the body spacethat is being mapped to the world space. R(t ) is the rotation matrix that gives the orientation. Conversionof rotation matrix to quaternion and vice-versa is discussed in Eberly (2010). Quaternions have a hugeadvantage over normal rotation matrices in a three dimensional space. Instead of storing 3×3 matrixindicating rotation using the rotation matrix, only four scalar values are needed if quaternions are used.The movement of the rigid body can be decomposed to linear velocity vc (t) of the center of mass andangular velocity ω(t) relative to the center of mass. ω(t) is the unitary axis of rotation multiplied by theangular velocity. The velocity v(t ) of a point that has world space position r (t ) is computed as

r (t ) = v(t ) = vc (t )+ω(t )× (r (t )− rc (t )) (2.5)

The linear momentum p(t ) is calculated from the linear velocity using the mass of the body

p(t ) = mvc (t ) (2.6)

and the angular momentum b(t ) is computed from the angular velocity and the inertia tensor

b(t ) = Iω(t ) (2.7)

Then during the time-step, the net force-torque pair acting on each of the objects are computed aswell as the other external forces acting on the objects. The rigid body’s movement is controlled by thelinear, and angular momentum applied. The change in linear momentum can be calculated as

10

∂p

∂t= F (2.8)

Where F is the force applied to the center of mass which changes the linear momentum of the rigid body.Similarly, to rotate the object torque, τ is needed. The torque is determined by the force F and the point ofits application r in the world space relative to the center of mass.

τ= (r − rc )×F (2.9)

the change in angular momentum is then given by

∂b

∂t= τ (2.10)

From the linear and angular momentums, the velocities are calculated and changes to the dynamicstate of the rigid body are updated. More detailed description can be found in Coutinho (2012).

During constrained motion of rigid bodies, the collision should be dealt properly. Otherwise, it will leadto unphysical results. So in each time-step, when multiple rigid bodies are involved, collision detectionis an essential step in analyzing the contact pairs as well as in the case of collisions resolving the forces.Modeling of collisions between rigid bodies and the analytical response has been discussed in Moore andWilhelms (1988). The discussion about collision detection between polyhedral rigid bodies by using theoriented bounding boxes (OBB) is discussed in Redon et al. (2002).

Testing for collision is often a costly process especially when the objects are irregularly shaped. Henceto reduce the complexity of this, the objects are normally approximated using a regular geometric primitive.Either it is approximated as a single primitive or union of the primitives using affine transformations.Common primitives that are considered are: spheres, cones, cylinders, boxes, and polygons. Theseprimitives are all convex shaped, meaning for every pair of points within the region, every point onthe straight line segment that joins the pair of points is also within the region. Convex objects oftenallow simpler and faster algorithms for intersection testing. More detailed description can be found invan den Bergen (1999). But not all bodies can be rendered as the union of particular primitive shapes toget a good approximation. Therefore developers normally use the discrete number of polygons whichmeans polygonal meshes to tessellate the rigid body surface. Most simple and commonly used polygon istriangles because of the following desirable properties.

• The triangle is the simplest of the polygon. Any fewer than three vertices and we would not have asurface at all

• A triangle is always planar because three points constitute a plane. More than three this may nothold true.

• Triangles remain triangles under most transformations including affine transforms.

• Almost most of the rendering platforms support triangle rasterization.

Because of the listed advantages, the crystal body used in the simulation has been triangulated using asurface of triangles. Rasterization techniques and triangulation are discussed more detailed in Gregory(2009). More of this will be discussed in section 3.1.

Checking for intersection in primitives is a time-consuming process, especially when it involves rigidbodies with a large number of polygons in the simulation. To minimize the computation costs objectbounding volumes are used. A bounding volume is a single simple volume encapsulating one or moreobjects or more complex nature. More on the bounding volumes and intersections tests can be found inEricson (2004). Once a collision is detected, it has to be handled accordingly. Two rigid bodies A and Bcollide at time tc . Let ra(tc ) and rb(tc ) be the point of contact of body A and B in world space concerningthe center of mass of A and B respectively. The points r A(tc ) and rB (tc ) coincide in the absolute coordinatesystem but their velocities (r )A(tc ) and (r )B (tc ) can be quite different. If the rigid body A is moving withlinear velocity v A and angular velocity ωA in time tc and analogously for B then with by equation 2.5

r A(tc ) = v A(tc )+ωA(tc )× r A(tc ) (2.11)

rB (tc ) = vB (tc )+ωB (tc )× rB (tc ) (2.12)

11

which gives the actual velocities of the contact points. The direction of the collision is described by theunit normal vector nB (tc ) of the surface of the rigid body B in point rB (tc ). This normal is pointing inthe outward direction of the volume of body B . To calculate the relative motion vr el of the two bodiesprojected in nB (tc ) direction

vr el = nB (tc ).(r A(tc )− rB (tc )) (2.13)

As shown in the Figure 2.1 based on the polarity of the relative velocity, the bodies may be inter-penetrating or moving away and if it is zero then they form a resting contact (Baraff (2001)).

Figure 2.1: Possible relative velocities of two objects (Kavan, 2003).

The collision point, collision normal and relative velocities are resolved for every contact detected inone time-step. Then collision-resolution module, which acts like a black-box in most rigid body simulationframework, deals with accordingly to prevent unphysical results. There are various approaches discussedby the researchers such as linear complementary problem approaches by Trinkle et al. (1997), penaltyapproaches by Donald and Pai (1990), Song et al. (2001) and, velocity-impulse LCP based time-steppingmethods by Stewart and Trinkle (1996) and Stewart (2000).

2.1.1 Physics Engine

Physics Engine pe is the highly parallelizable multi-body simulation framework originally developed byKlaus Iglberger (Iglberger and Rüde, 2009) at the University of Erlangen-Nuremberg. Geometric primitivesare utilized so that complex shapes can also be used in the simulations. The rigid body type developed inthis work is also works as an example of how a complex shape can be implemented in pe.

pe uses Message Passing Interface (MPI) parallelizations to carry out large-scale simulations. Thedomain is divided into multiple sub-domains, and each sub-domain is operated upon by a separateprocess. All the process knows their neighboring process, i.e., when they share a common boundary witheach other, and at each time-step, they communicate with each other, about bodies intersecting in theconnected domains and ownership of moving bodies. Communication and ownership modules used inpe can be found in Iglberger and Rüde (2009) and in Preclik (2014)

Impulse/force approach had been implemented initially in pe to solve the collision problem (Iglberger,2010). This method involves merging all the contacts into a single global problem. Then a Fast FrictionalDynamics (FFD) algorithm was implemented to improve the scalability of the solver and that was improvedby Wengenroth (2007). In this solver, each collision is treated locally. Sometimes FFD algorithm mightnot be physically accurate (Götz, 2012) due to the simplification of some contact scenarios. To overcomethis oversimplification, a time-stepping solver for hard contacts called the Hard Contact Semi-ImplicitTime-Stepping solver (HCSITS) was implemented in pe by Preclik (2014). For the scope of this work, thissolver is used.

Using an efficient and effective method to find the intersection between two bodies is imperative. Onesuch algorithm is the Gilbert-Johnson-Keerthi Algorithm developed by Gilbert et al. (1988) to find theintersection between two convex polyhedra. It uses a simplex-based method to find the closest points ofthe convex set of points. The separation distance between the two convex objects is calculated using theMinkowski sum or also known as Minkowski difference. More on this algorithm and Minkowski distancecan be found in Ericson (2004).

The basic algorithm involving pe is given as follows (Preclik, 2014)

12

Algorithm 1 pe Rigid body simulation time step

1: procedure SimulateTimeStep(δt )Require: All Shadow Copies and Global Bodies are in syncRequire: Shadow Copies must be available on all intersecting subdomainsRequire: All bodies must be located in their owners’ subdomains

2: Detect Contacts3: Filter Contacts4: Initialize velocity accumulators(δt )5: Synchronize velocities6: Contact Reaction (Collision Response)7: Synchronize velocities8: Integerate Positions(δt )9: Sync Positions

Ensure: All Shadow Copies and Global Bodies are in syncEnsure: Shadow Copies must be available on all intersecting subdomainsEnsure: All bodies must be located in their owners’ subdomains

10: end procedure

Collision detection is followed by determining an appropriate response. For example, when a stone isthrown against the wall it should hit the wall and bounce off when it is assumed both the bodies are rigidand cannot undergo deformation. At the point of collision, the bodies change direction and velocity.

Contact generation plays an important part to correct any and all interpenetration that rigid bodiesundergo. The following parameters describe all the detected contacts.

• Rigid bodies in contactThe two bodies can be referred to as body A and body B for the scope of this explanation. The

bodies will either are in contact with each other or inter-penetrate each other to initiate a collisionresponse.

• The point of contactThe intersection between the two points will generate a point which is known as a point of

contact. This point of contact will be represented in world frame of coordinates. This point cannotbe assumed to be on the surface of one body instead it can be any point in the intersection. Thispoint of contact will be generated during the detecting contacts step and will be stored in a containerfor further action.

• The contact normalThe contact normal conventionally points from the center of mass of body B to body A’s center

of mass. Usually a unit vector, contact normal affects the algorithm and the mathematics used toresolve the contact point by applying suitable penetration correction.

• The contact distanceThis distance indicates the gap between the two bodies undergoing contact or about to undergo

collision. If it is a positive value, then there is no inter-penetration, or there might not even be anycontact. If it is a negative value, then the bodies have inter-penetrated. The contact distance gives avalue that is the negative of the penetration distance. This penetration depth can be used by thecollision response algorithm to calculate the appropriate response in the direction of the contactnormal.

The Figure 2.2 shows the various collision parameters that were described above. It shows the penetra-tion distance, contact normal and contact point between two bodies.

13

Figure 2.2: Collision Parameters (Schornbaum, 2010).

Once the penetration is observed, the suitable penetration correction is prompted. This correction isonly executed when the penetration depth is lesser than the “Penetration Threshold”. If it is more, thenthe penetration correction does not affect this contact. To take into account the instability and inaccuracycaused by floating point arithmetic, the penetration threshold is set as close to zero as possible andtherefore have a highly accurate contact detection. Once the contact and penetration are observed, it hasto be resolved. The idea behind it is to apply proper impulse based on the linear and angular velocities ofthe bodies before the collision. This procedure has to conserve the velocities before and after the collision.The factors affecting the collision response are the type of collision, elastic or inelastic and also the frictionbetween the materials which is controlled by the coefficient of friction. The effect of these parameters onthe collision response is discussed in Schornbaum (2010) and Preclik (2014).

2.2 Lattice Boltzmann MethodIn the field of computational fluid dynamics, solving the problem numerically involves solving the

continuity and the Navier Stokes Equation (NSE). In this approach, the fluid is considered as a continuummedia leading to the fluid to be considered a continuous blob of matter rather than a collection ofmolecules. This scale is known as macro-scale. If one considers only the inter-molecular forces of the fluidthat were interacting with each other then according to Newton’s second law of momentum, this leads tothe Molecular Dynamics simulation. This scale is known as micro-scale. In between these two scales existthe mesoscale approach, which deals with the distribution of particles. This mesoscopic kinetic theoryforms the crucial component in Lattice Boltzmann Method (Krüger et al., 2017).

The Boltzmann equation can be given as

∂ f

∂t+ξβ

∂ f

∂xβ+ Fβρ

∂ f

∂ξβ=Ω( f ) (2.14)

where f represents the distribution function f (x,ξ, t ) , ξβ the velocity vector of particles . Ω represents thecollision term which influences the distribution function based on the collisions. The original collisionoperator is difficult and rather cumbersome because of the double integral over the velocity space. Thesimple and the most commonly used discretization method for the collision operator with single relaxationtime was developed by Bhatnagar et al. (1954). The so called BGK collision operator developed byBhatnagar, Gross and Krook is as follows

Ω( f ) =−1

τ( f − f eq ), (2.15)

where τ is called the relaxation time. This relaxation time captures speed of the relaxation of the distribu-tion function towards the equilibrium distribution( f eq ). Discretization of the Boltzmann equation leadsto the discrete Lattice Boltzmann model (Qian et al., 1992). By choosing q number of discrete velocities ci ,the discrete LB equation can be written as

∂ fi (x, t )

∂t+ ci .∇ fi (x, t ) =Ω( fi ) (2.16)

and the macroscopic moments density ρ and momentum are also computed as follows

14

ρ =∑i

fi =∑

if eq

i (2.17)

ρu =∑i

fi ci =∑

if eq

i ci (2.18)

and pressure is related to density and can be obtained by a simple relation as follows

p = c2s ρ (2.19)

where cs is the speed of sound. The first order discretization of LB equation 2.16 in space and time usingthe finite differences leads to

fi (x + ci∆t , t +∆t )− fi (x, t ) =∆tΩi (x, t ) (2.20)

where ∆t is the time-step. Substituting the equation 2.20 in 2.15 gives us the Lattice Boltzmann Equationwith BGK collision operator

fi (x + ci∆t , t +∆t ) = fi (x, t )− ∆t

τ( fi (x, t )− f eq

i (x, t )) (2.21)

this equation 2.21 is also called as LBGK equation. From this equation one can separate two differentsteps, namely collision step and streaming step.

f ∗i (x, t ) = fi (x, t )− ∆t

τ[ fi (x, t )− f eq

i (x, t )] (2.22)

fi (x + ci∆t , t +∆t ) = f ∗i (x, t ) (2.23)

where f ∗i is called the post collision population. The collision step is responsible for each population to get

the collision contribution and edge towards the local relaxation, whereas the streaming step is responsiblefor the post collision population to the neighboring lattice cells.

Discretization of velocities of the LBM equation will lead to choosing of specific lattice model becausechoosing the right lattice model is vital so that the NS equation is well resolved and it is not computationallyexpensive. The lattice model is denoted as DdQq where d is the dimension, and q is the number of discretevelocity sets. D3Q19 has been chosen because of the compromise in efficiency and accuracy (Mei et al.,2000). This lattice model is shown in Figure 2.3 which has 19 PDFs, and the 19 discrete velocities of theneighboring lattice cells can be formulated as (Feichtinger, 2012)

Figure 2.3: D3Q19 Model (Iglberger et al., 2008).

15

ci =

(0,0,0) i ∈C(±c,0,0) i ∈ E ,W(0,±c,0) i ∈ N ,S(0,0,±c) i ∈ T,B(±c,±c,0) i ∈ N E , NW,SE ,SW(±c,0,±c) i ∈ T E ,T W,BE ,BW(0,±c,±c) i ∈ T N ,T S,B N ,BS

(2.24)

and by applying Taylor expansion to the local Maxwell distribution f eq in equation 2.21 for low Machnumbers, for isothermal flows and for incompressible fluids (ρ0 = 1) leads to

f eqi (ρ(x, t ),ui (x, t )) = wi

[ρ+ρ0

(ci u

c2s

+ (ci u)2

2c4s

− u2

2c2s

)](2.25)

where the lattice weights wα for the D3Q19 model can be given by (Succi, 2001)

wα =

13 if ci = (0,0,0)

118 if ci = (±c,0,0), (0,±c,0), (0,0,±c)

136 if ci = (±c,±c,0), (±c,0,±c), (0,±c,±c)

(2.26)

where the local equilibrium function is dependent on macroscopic density ρ(x, t) and macroscopicvelocity u(x, t ) which can be calculated using equation 2.17. Also from equation 2.19, pressure for D3Q19can be simplified to

p = 1

3ρ (2.27)

because relation between lattice velocity c and the speed of sound cs is given by

c = ∆x

∆t= 1 (2.28)

cs = 1p3

c (2.29)

Based on the relaxation time τ various collision model variants exist. Previously we have discussedSingle-Relaxation time (SRT) variant, but it has stability problems, especially when viscosity is small andalso when no slip boundary conditions are used. When using no slip boundaries, the slip velocities at thewall depends on the relaxation parameter τ. This results in the virtual displacement of the boundaries Götz(2012). In order to overcome these stability issues, Two-Relaxation Time (TRT) that has been developed byGinzburg et al. (2010) is used. This method is comparable to the SRT method in terms of complexity butmore stable than the SRT.

2.2.1 waLBerla

Widely applicable Lattice-Boltzmann from Erlangen (waLBerla) is a specifically developed frameworkat the Chair of System Simulation at the University of Erlangen-Nuremberg to carry out massively parallelfluid flow simulations centered around the Lattice-Boltzmann method. More detailed description ofwaLBerla can be found in Feichtinger et al. (2011) and in Feichtinger (2012).

The implementation of waLBerla includes distributed problem setup meaning that the domain willbe split into blocks and each block will be operated upon separately. The data allocated in blocks isindependent of the other. In waLBerla all the operation that is carried out during one time-step is alsosplit into parts known as the sweep. During each sweep, the data in the grid are worked upon. Even thedata fetching and data writing are split into one of these sweep routines.

The LBM sweep applied on the simulation domain can be given as

16

Algorithm 2 LBM Sweep

1: procedure SimulateTimeStep(δt )2: for <every block> do3: Setup Boundaries // setup boundary conditions4: Stream step //propagate information through the grid5: for <every cell> do6: Calculate macroscopic variables // e.g : density7: Calculate the equilibrium distribution function from macroscopic varibales8: end for9: Collide step // relaxation towards function distribution function

10: end for11: MPI Communication // communication between blocks12: end procedure

2.3 Coupling - Momentum Exchange MethodLattice Boltzmann method is highly efficient algorithm to solve fluid flow problems. Accordingly

when particles are also included in the simulation, the method used to couple the particle flow shouldalso be numerically efficient. To have efficient fluid-solid coupling to carry out the particle-laden flows,momentum exchange method is being used. It was first proposed in Ladd (1994).

The idea behind this method is to map cells containing the solid body in the flow field explicitly asobstacle cells and to apply the no-slip conditions at the cells on the surface of the solid. No-slip boundarycondition, because to make sure that the fluid around the surface of the particle is moving at the samevelocity as the particle. Ladd (1994) used the bounce-back boundary condition to compute the fluidforce on the particle. The particle was discretized into some cells, effectively turning it into a staircaseapproximation. The mapping of the obstacles and algorithm for the treatment of particles in the fluid flowsimulation will be discussed in the next section 6.

2.3.1 pe - waLBerla Coupling

In this section, simulation method that is used to couple the rigid bodies with fluid flow is addressed.The LBM method used in waLBerla is extended by coupling it with pe to resolve the motion of particlesand their interaction with the new rigid body type.

The coupling of the two frameworks has to take care of a couple of issues. The first one is the mappingof the rigid bodies to the flow simulation. Mapping can depend on the rigid body on the basis that it ismoving or not. For a moving rigid body, the boundaries have to be mapped continuously. The next issue isthat to compute the hydrodynamic force acting on the body.

The particle treatment that is happening in each step is mentioned below (Götz, 2012)

• Mapping of the obstacle cells and reconstructing the PDFs when the particle moves

• LBM stream and collide step

• Force Evaluation

• Movement and collision of objects

Moving objects in the flow simulation have to be mapped in both these frameworks in different ways.In pe the rigid bodies are stored as a Lagrangian representation. In waLBerla the objects are mapped in aEulerian frame, meaning the objects are mapped to the lattice grids. The different representation of theobject is shown in Figure 2.4.

In pe, the sphere is a complete body with a center of mass and a radius. In waLBerla the sphere ismapped to the grid by computing the cells which have the center inside the sphere and assigning themwith the obstacle flags. This leads to discretization error as shown in the Figure 2.5.

When the object moves during the simulation at any time-step, from its previously mapped area in thegrid, the obstacle has to be remapped, i.e., the cells that were previously fluid has to mapped as obstaclesand the cells that the object has moved away from has to be reconstructed with the missing PDFs whichhas been explained in Pickl et al. (2011), Götz (2012) and Rettinger and Rüde (2017). The next step is theforce evaluation, i.e., calculating the hydrodynamic forces that are acting on the rigid body during the

17

Figure 2.4: Representation of Objects in pe and waLBerla. On top : Simulation setup with sphere in achannel flow . Lower left : representation in the LBM fluid solver. Lower right : representation in the rigidbody solver (Götz, 2012).

Figure 2.5: Left : Discretization error for a sphere mapped to LBM grid . Right : Zoomed in figure of thediscretization error (Götz, 2012).

streaming step. When the fluid particles stream from their cells, the cells occupied by the solid (obstaclecells) will revert the process resulting in the momentum exchange between particles and fluid (Pickl et al.,2011).

Algorithm 3 shows the algorithm used for the coupling (Rettinger and Rüde, 2017).

Algorithm 3 LBM Algorithm with Momentum Exchange Method

Initially map bodies to the fluid domain1: for <each time-step t> do2: for <each LBM sub-cycle> do3: Perform LBM collision step4: Apply the boundary conditions5: Perform LBM stream step6: Calculate the hydrodynamic forces on particles7: end for8: Average the forces on particles over LBM subscycles9: Obtain total force and torque on particles

10: for <rigid body solver subcycle> do11: Perform rigid body solver step (collision detection and resolution and time integeration)12: end for13: Update particle mapping and reconstruct PDFs if necessary14: end for

18

3 Domain setupIn this section we will discuss how the creation of the rigid body is carried out and how it is imported

in pe and in waLBerla.

3.1 Crystal Body3.1.1 Convex Polyhedra

Geometric primitives are the building block of the simulation environment. However, when theshapes become complex, then the approximation using the primitives becomes difficult. So in the fieldof computer graphics use of polygonal meshes have been prominent. Polygonal meshes are surfacerepresentations, which reduces three-dimensional problems into two-dimensional problems.

A polygon is a closed figure with n-sides, and each side is a representation of connected points in aplane. Each point in the polygon is called a vertex and the line connecting the vertex is called an edge. Apolygon is called convex when the line connecting any two points inside the polygon remains inside thepolygon. This is shown in Figure 3.1

(a) Polygon and its geometrical segments (b) Convex Polygon

Figure 3.1: Convex Polygon (Ericson, 2004).

The complexity of the triangle meshes is measured by the number of vertices and number of faces. Atriangle mesh M can be represented by a set of vertices

V= v1...., vV (3.1)

and a set of triangular faces connecting them

F= f1....., fF, fn ∈V×V×V (3.2)

where each triangle specifies its three vertices from V. Each vertex is associated with a 3D position pi in R3

and can be specified as

P= p1.....,pV, pi = p(vi) =x(vi)

y(vi)z(vi)

∈R3 (3.3)

such that each f ∈F represents a triangle specified by its three vertex positions (Botsch et al., 2006). An-other important parameter to describe the triangles in the mesh are the edges. This gives the connectivityof the vertices of the triangles. It is given by

E= e1.....,eE, ei ∈V×V (3.4)

So consider a triangle with points p1,p2, and p3 which denote the vertices. Then from these points,the edges can be computed using

e12 = p2 −p1, (3.5)

e13 = p3 −p1, (3.6)

e23 = p3 −p2. (3.7)

19

and the face normal of the triangle can be calculated by

N = e12 ×e13

|e12 ×e13|(3.8)

and it can be represented as shown in Figure 3.2.

Figure 3.2: Edges and Face Normal of a Triangle Gregory (2009).

In this work, only the non-deformable bodies are considered. To describe a non-deformable object inthe physics-based simulation, only the surface is approximated with a polygonal mesh. This procedure ofapproximating using a surface with a piecewise linear representation is known as tessellation. Using a2-simplex triangular mesh for this tessellation procedure is known as triangulation Gregory (2009).

These triangular meshes are used extensively to approximate the convex and non-convex polyhedrons.A polyhedron is three dimensional body formed by connecting a set of plane polygons to form a closedsurface, which splits the polyhedron into two regions, one inside the polyhedron called the interior regionand other outside the polyhedron called as the exterior region. The convexity of the polyhedron canagain be characterized similarly to the convexity of the polygon. A polyhedron is convex if all the pointsin the line joining the two points inside the polyhedron is contained in the interior region. A convexpolyhedron can be described as the intersection of a finite number of half-spaces. More on the half-spacesand simplices can be found in Ericson (2004). This description of interior and the exterior region will helplater in detecting collisions with spheres.

3.1.2 Creation and Import

The crystal structure design used in the simulation is constructed in such a way that it imitates thestructure of the diesel particulate filter made from the advanced ceramic material based on acicularmullite. This filter structure consists of long needle-like structures randomly distributed as shown inFigure 3.3 which creates a highly porous diesel filter, and they have extreme mechanical integrity. Highporosity also results in the ability of the filter to retain a lot of catalytic materials but also not compromisingthe back pressure. More on the characteristics of the ACM filter and its material properties can be found inPyzik and Li (2005).

Figure 3.3: Texture of ACM Diesel Particulate Filter Pyzik and Li (2005).

20

From the Figure 3.3, it can be seen that these needle-like components have somewhat unevenlysized thickness and length. Moreover, the position of the filter, as well as the orientation, is also random.Parameters that controls the creation of the crystals are shown in the Figure 3.4. Actual thickness andlength scales used in the simulation vary based on the type of simulation.

Figure 3.4: Single crystal structure used in the filter.

Creation of the crystals is assumed to be from the material whose dimensions are known. Crystals aresaid to be developed from a point known as the seed located within the crystal and do not grow the samein both directions, resulting in the eccentricity which is shown in the Figure 3.5.

Figure 3.5: Eccentricity of the Crystal (Gil et al., 2017).

After the creation, rotations are applied to the elevation and azimuth creating an uneven and chaoticcrystal structure as shown in the Figure 3.6.

Figure 3.6: Base Crystal Structure.

21

Crystals are defined using 16 triangles, that encapsulates 10 points and 24 edges (Gil et al., 2017). Inthe Figure 3.4, it is shown how one side is split into four triangles, and it is repeated on the other five sides.Then the crystals are exported as STL(STereoLithography) file. STL files are widely supported in varioussoftware platforms and are prominently used to represent surface or body, that is tessellated using trianglemeshes. It contains the triangle vertices of all the crystals along with the face normal. All the vertex pointsof the triangles have a common origin and describe the position of the crystal in the simulation domain.

The next step is to read the data from the STL file and to create the rigid body type by the name ofCrystalBody. In pe all the rigid body types are represented under the class GeomPrimitve similarly to theother geometric primitives like a sphere, box, and cone. Then the new rigid body type implemented canuse functions such as setPosition, setCommunicating and registerContact which is defined in the mainclass.

Before discussing about how crystals are created as rigid bodies, how the data is imported from the STLfile has to be discussed. For this, the class hierarchy and the objects created has to be examined. Duringthe parsing of the file, three classes are essential which are already implemented in waLBerla. First is theTriangleMesh class, which can be used to read and write triangle meshes. However, the two main classesare CrystalSet and Crystal. The dependency is shown in the Figure 3.7.

CrystalSet class is called from the actual source file with the STL file being loaded from the parameterfile. This class, while reading the file, will distinguish between each crystal separately. The crystals canbe translated to any position in the domain with the help of the routine implemented in this class. It iscalled everytime a crystal is created, because all crystals are moved to the center of the domain ensuringthe formation of a wall of filter crystals.

Figure 3.7: Class Dependency.

Once a single crystal is identified from the mesh file, the Crystal class object is created. This classobject will contain the details of all the triangles that describes this crystal. Scaling of the crystals is alsodone in this step while creating the crystal object. The mesh is used to map the boundary of each crystal byforming an Axis Aligned Bounding Box (AABB). This bounding box method is one of the common methodsto describe the bounding volume of the rigid bodies. AABBs are used in pe framework Iglberger (2010) toaccelerate the collision detection routines. The TriangleMesh is called from the constructor of the Crystalclass which stores the vertices of all the triangles in an array structure.

Till now the process of reading the STL mesh file and creating the crystal object was described. Thenext step takes this crystal object and creates a rigid body in the pe framework. In pe, the creation of rigidbodies involves specifying the material and also some rigid body specific properties that are used later inthe simulation and they are explained as follows

• Global: A boolean value, specifying whether the given body is global or not. Switching on the globalparameter will set the flag so that this rigid body is visible to all the processes in an MPI simulation.

• Communicating: Setting this flag will make sure that the neighboring process will know the presenceof the object. If the rigid body is sliding into another block, one who does not hold the ownership of

22

the said rigid body, this flag will make sure that a shadow copy of the rigid body is created in theneighboring process.

• InifinteMass: Setting this flag will make the rigid body fixed.

More about the shadow owners and the communicating procedures used in pe can be found in Iglberger(2010). For this work the flags used in the simulation are as follows

GlobalThis flag is set to false. The reason being each crystal is a process local body. During the creation of

the crystal body, the centroid of the crystal is calculated. The centroid calculation involves calculatingthe barycenter of all the triangles. Let xi , yi , zi be the vertices of triangles. Then the barycenter, Bi can becalculated as

Bi = xi + yi + zi

3(3.9)

Then the area of the triangle, Ai is calculated using

Ai =∥∥(yi −xi )× (zi −xi )

∥∥ (3.10)

Then the centroid, Ci can be given as (Nürnberg, 2013)

Ci =

N−1∑i=0

Ai Ri

N−1∑i=0

Ai

(3.11)

Once the centroid is calculated, the process or the block as it is called in pe checks whether the centroid iscontained in it. If the centroid is present in that block, then the crystal body is created in that block, andthe ownership belongs to that process.

Infinite MassThis flag is set to true. Setting infinite mass to true makes sure that the simulation cannot move the

rigid body during the position integration of rigid bodies. The crystals are supposed to be fixed at thecenter of the domain.

CommunicationThis flag is set to true. The crystals are usually spanned more than one block. Therefore the neighbor

process should have a shadow copy of the crystal. The idea of the shadow copies and parallelizationtechniques are available in Preclik (2014) and Iglberger (2010).

Along with the simulation parameters required to describe the rigid body the material was also chosen.The material can be stipulated by two parameters namely static and dynamic coefficient of friction. Theseparameters play an essential role in the collision response observed from the rigid body.

Creation of the rigid body is also followed by the process of syncing and communication of all theprocess, so that information of all the bodies in the neighboring processes are appropriately disseminated.Once the creation is done, the actual mesh file containing all the crystals is written as an OBJ file. This willhelp us in visualizing the actual crystal filter structure used in the simulation. A prototype of the samplefilter is shown in the Figure 3.8.

3.1.3 Collision Detection between Crystal Body and Sphere

The next step in the setting up the new rigid body is the process of setting up the collision detectionand response. During the simulation, the spheres might not be near the crystal body to undergo collisionin the next time step. This process of elimination based on the impossibility of collision to happen isknown as broad phase collision detection. Broad phase collision detection can be implemented usingspatial partitioning.

In spatial partioning, the domain is split into cubic domains and all the rigid bodies, which is encapsu-lated by the bounding volumes, are assigned to the grids over which it spans. Rigid bodies that are not inthe same grid or neighbouring grids are not tested for collision as they have no possibility of undergoing

23

Figure 3.8: Prototype - Filter Simulation.

collision in the next timestep. More on the use of this method to improve coarse collision detection can befound in Schornbaum (2009) and Iglberger (2010).

The next step is to check for collisions between bodies that were filtered from the broad phase collisiondetection using Bounding Volumes. They are easy to construct with minimal memory requirementsand could contain the rigid body in their entirety, meaning that it will cover all the extremes in all threedirections. In this work, AABBs are used. During the narrow phase testing, the bounding boxes for therigid bodies are checked against each other for an intersection. If there is an intersection, then a possiblecontact pair is generated, which will be dealt with later. It accelerates the chance to rule out a potentialcollision between pairs. More about the bounding volume hierarchies and also about the intersectiontests can be found in Ericson (2004).

From these possible contacts, the actual contact between pairs is checked. This continuation of narrowphase collision detection depends hugely on the primitives that are being tested. Usually for convexobjects, the Gilbert-Johnson-Keerthi (GJK) algorithm is used (Ericson, 2004). However, the implementationused here determines the penetration using analytical methods genertating the contact parameters. Thecollision between the two primitives is checked based on the Collide function that is implemented in theRigid body class, here CrystalBody. In this work, only the routine to check for the collision between theCrystalBody and sphere is implemented.

Figure 3.9: Voronoi Region - 2D triangle (Scharpff, 2014).

So the actual implementation involves testing the sphere against all the triangles of the rigid body andgenerating the closest point between the sphere and the CrystalBody. The first step is the get the positionand the radius of the sphere. Then the position of the sphere is checked against each triangle to find theclosest point. While checking for the closest point, the triangle is split into eight regions known as Voronoiregions as shown in Figure 3.9.

The algorithm 4 to check for the closest point to the triangle involves finding the Voronoi region inwhich the point lies and then projecting the point towards that Voronoi region. The implementation isalready found in Ericson (2004) and also given in Scharpff (2014).

24

Algorithm 4 Finding the Closest Point to Triangle

1: procedure sqDistance(Point O, Point ClosestPoint, Region)Vertices A, B, C

2: Compute−→AB ,

−→AC ,

−→AP

3: if (−→AB ×−−−→

ABC )T ×−−→AO > 0 then

4: if−→AB T −−→AO > 0 then

5: return region = 2

6: else if−→AC T −−→AO > 0 then

7: return region = 58: else9: return region = 1

10: end if11: else12: if (

−−−→ABC ×−→

AC )T −−→AO > 0 then13: go to Row 6 . Same as this code to get regions 7 and 8 respectively

14: else if−−−→ABC T −−→AO > 0 then

15: return region = 316: else17: return region = 418: end if19: end if20: end procedure

Once the region is found out, then the Closest point is returned which is stored for all the trianglesof the crystals that are shortlisted from the narrow phase collision detection. Then it will be checkedfor penetration, meaning that when the closest distance is less than the contact threshold, then the twobodies are in contact. The algorithm 5 for this function is as follows

Algorithm 5 Collide Function

1: procedure Collide()2: for All crystal Bodies from narrow phase collision detection do3: for All triangles do4: Closest Distance() . Closest point and the normal from Algorithm 45: end for6: if Closest Distance < (Contact Threshold + Sphere Radius) then7: generate contact and contact normal and add it to the container8: end if9: end for

10: end procedure

The crystal body checks all the triangles for a possible collision. If there is a collision happening, itfinds out the penetration distance and returns that as the contact distance between the sphere and thecrystal body. At the same time, the contact normal is also constructed. Then an impulse is triggered basedon the linear and angular velocity of the sphere and applied in the direction of the normal, in this casepointing from the crystal body towards the sphere. Since the crystal body is fixed the impulse force acts onthe sphere.

3.1.4 Mapping in waLBerla

Coupling of the LBM solver and rigid body dynamics simulation involves mapping the rigid bodies tothe lattice cells. If the midpoints of the lattice cells are inside the rigid bodies, then they are considered asobstacle cells or moving boundaries based on which body is considered. Since the filter is a fixed body, itdoes not need to be mapped after every time-step as the sphere. The primary criteria is whether the centerof the lattice cell is inside the crystal or not. To find whether a point is inside the triangular mesh, themethod of ray intersection test is used. The idea is to create a ray from the centroid of the crystal structure

25

towards the point and when the intersection query is generated, the point is inside the structure if thevector pointing from the centroid to this cell center did not cross the triangle. If the point is outside thetriangle then the vector pointing from the centroid will cross one of the triangle faces or a vertex.

Figure 3.10: Mapping of Crystal to waLBerla.

The crystal structure to be mapped is checked against all the cells in the domain. Then the vectoror the ray pointing from the centroid to the cell to be checked is constructed. This vector is checked forintersection against all the triangles that describe the crystal. The first check is whether the vector istraveling towards the face, by checking against the face normal of the triangle. The check shown in step 11will find out whether the vector pointing from the centroid towards the cell center is in the same directionas the face normal. If it is either parallel or in the opposite direction, it will exit.

The next checks will determine whether the vector from the centroid intersects with the face of thetriangle. Checks in step 16 and 18 will determine the intersection of the face and check 24 and 26 willmake sure the barycentric coordinates are within the bounds. The algorithm can be found in Ericson(2004) and also a detailed explanation of projection of vectors and barycentric coordinates is also given. Ifall the tests are cleared, then the vector to the Cell center from the centroid is intersecting, meaning thepoint is outside the crystal body. Then this cell can be marked as a fluid cell. This routine is repeated forall cells against all crystals, and the cells which are inside the crystal body are mapped as no-slip cells.

26

Algorithm 6 Mapping Of Crystal Body to waLBerla

1: procedure CrystalBody::containsRelPointImpl()2: for All the lattice cells do3: P ←center point of the lattice cell4: for All the crystal Body do5: C ←Centroid of the CrystalBody6: for All the triangles do7: A,B ,C ← vertices of triangle

8: Compute−→AB ,

−→AC and

−→C P

9: compute face normal, ~N =−→AB ×−→

AC

10: D =−→C P T N

11: if D < 0 then12: return false and exit . False - intersecting the face - Point is outside13: end if14: compute

−→AP = P − A

15: T =−→AP T N

16: if T < 0 then17: return false and exit18: else if T > D then19: return false and exit20: end if21: E =

−→C P ×−→

AP22: V =

−→AC T E

23: W = -−→AB T E

24: if V < 0 || V > D then25: return false and exit26: else if W < 0 || V+W > D then27: return false and exit28: end if29: end for30: return true . If all tests fail, Not intersecting - so Inside the Crystal body31: end for32: end for33: end procedure

27

4 Sizing

4.1 Crystal Body SizingWhile importing the mesh file to create crystals as rigid bodies, a new problem is discovered. The

length of scale of crystals is in micrometer. The minimum length of the crystal is 10µm and the maximumis 100µm. The thickness to length ratio ranged from 0.01 to 0.05. This means the thickness ranged from0.1µm to 5µm. The crystal mesh file that was generated based on these parameters did not follow theactual physical scale.

The domain used for simulations in pe are adapted to match the length scales of the crystals from themesh file. Therefore, there is no need to resize the rigid body. Particles are described using radius andthe position alone. But in waLBerla the default size for the spatial discretization in the lattice unit is setto 1 (∆x∗ = 1). So when transforming from the physical units to the lattice units, the domain includingthe rigid body has to be scaled. This means when the mesh file is imported in the simulation, the entirecrystal set has to be scaled to so that it can be mapped. Scaling process also influences the domain size.

Let Lpmax be the maximum length of the crystal in some physical unit in the mesh file. The actual

maximum length of the crystal found in the crystal set is 100µm. Let this be Lamax . let Lp

mi n and Lami n

minimum length in mesh file and in actual length scale respectively. Lami n is already known, 10µm.

When measured the maximum and minimum is found out to be

Lpmax = 0.103762 (4.1)

Lpmi n = 0.0108315 (4.2)

from this value the scale factor is calculated by

Scal e F actor = Lamax

Lpmax ∗∆x

= 100∗e−6

0.103762∗1e−6 = 963.74 (4.3)

from this scaling factor, the size of the base crystal structure 3.6 is found out as

Leng th = 196.29µm (4.4)

Br ead th = 104.989µm (4.5)

Depth = 134.25µm (4.6)

and the lattice cells it is mapped to be is also calculated as if ∆x∗ is set to 1

Leng th = 196 cel l s (4.7)

Br ead th = 105 cel l s (4.8)

Depth = 134 cel l s (4.9)

where ∆x is the space discretization assumed to be taken, and the ∆x∗ is the lattice unit cell size. Thescaling of the crystals is done by manipulating the vertices of the triangle. The scaling of the mesh wasalready implemented in waLBerla. When the mesh file is read, and the crystal class object for all thecrystals are created and before they are imported into pe the scaling process is taken care of. While passingthe configuration file or the input file, the scaling parameters are passed on.

4.2 Particle Sizing and Boundary ConditionsLBM simulations are mostly performed on lattice units. Lattice units are the non-dimensionalisation

of the physical parameters using dimensionless numbers. In waLBerla, the lattice units conversion mustbe done pre-simulation manually and passed on as simulation parameters. The transformation of physicalparameters into lattice units is an important step in maintaining stability, accuracy, and efficiency.

The physical domain can be converted to the lattice units by a conversion factors. The law of similaritystates that two systems are the same if they have the same Reynolds number Krüger et al. (2017). Particlescan be classified into two categories namely spheres and agglomerates. The spheres are typically 10−2µm.

28

The agglomerates are a collection of these spheres and range from 10−1µm to 1µm. But for the scope ofthis work only spherical particles are considered. If particles considered are 10−2µm, then scaling factorshould be calculated appropriately.

Fixing the number of cells describing particles to ten, the resolution of ∆x is 10−8µm. The scalingfactor is recalculated, giving us

Scal eF actor = Lamax

Lpmax ∗∆x

= 100∗e−6

0.103762∗1e−8 = 96374.39. (4.10)

which will lead to a size of the base crystal structure mapped to the lattice units can be calculated as

Leng th = 19644 cel l s (4.11)

Br ead th = 10506 cel l s (4.12)

Depth = 13435 cel l s (4.13)

Including the fluid domain in the inlet and outlet, leads to an unreasonable simulation size whichrequires a lot of computation time. Hence, the ∆x is set to 1e−6. This means the diameter of particles isalso increased to 10µm. Next step is to set the length of the domain and from that calculate other latticeparameters needed for the simulation. The depth of the filter is taken as the basis from which the domainis expanded. The filter is placed at the center of the domain, and then the domain is expanded on eitherside with the same measurement. The domain is showed in the Figure 4.1.

Figure 4.1: Domain size in waLBerla.

The simulation parameters such as lattice inflow velocity and the lattice viscosity are calculated bysetting the Reynolds number and relaxation time. The Reynolds number is chosen as 15 and thereforefixing the inflow velocity which in turn fixes the lattice flow velocity. The relaxation parameter τ∗ is set to0.9091, which is used to calculate the lattice viscosity ν∗ as follows

ν∗ = (c∗s )2(τ∗ − 0.5) = 0.13619 (4.14)

where c∗s is the lattice speed of sound which is given byp

1/3 (≈ 0.577). Then the inflow velocity can becalculated as

U = Re ∗νD

= 15 ∗ 6.1e−7

400∗1e−6 = 0.02287 m/s (4.15)

U∗ = Re∗ ∗ν∗D∗ = 15 ∗ 0.13619

400= 0.00511 (4.16)

where U∗ is the lattice inflow velocity and U is the physical velocity.

29

Figure 4.2: Domain with Filter

The boundary conditions set in the domain influences the flow of particles. The inflow boundarycondition is set so that the fluid is driven with constant velocity. So on the other face, the condition is setas pressure outflow. This means the fluid will be at atmospheric pressure. Then the next thing to consideris the walls. To have a simplistic model, no-slip boundary conditions are used. However, to prevent thespherical particles going out of the domain rigid body plane along the borders of the domain were created.They behave like the wall, meaning that particles can bounce off of it.

30

5 Simulation ResultsIn this section, the simulation that was carried out and their results will be discussed. First, the

interaction between the spherical particles and the crystal body is alone studied without the influence ofthe fluid. This is done in pe alone. This step will emphasize on the interaction between the sphere and thecrystals and various checks to make sure that they are working perfectly fine. All the simulation are carriedout using Hard Contact Semi-Implicit Time Stepping (HCSITS) solver implemented by Preclik (2014).

Then the next step is analyzing the results after importing the crystals inside the waLBerla framework.The fluid drives particles through the filter crystals, and some get trapped. In this part, the influence of thefilter structure especially the porosity of it is studied. Also, the relation between the particle size and thecrystal size is also studied. Although particles are not realistically sized, this study will still indicate howparticles get trapped in the filter.

All the simulations are carried out on the high-performance cluster maintained by the Department ofSystem Simulations (LSS) at the Friedrich-Alexander-University Erlangen-Nuremberg. The cluster consistsof eight compute nodes, with each compute node consisting of 32 cores. The simulations are carried outbased on the problem size on multiple nodes or on a single node.

5.1 Results with peBefore performing the simulations with particles in pe, the rigid body primitive that was implemented

has to be tested. For this various simulations are carried out, and the results are analyzed. After theimplementation, visual check of spheres interacting with the crystal body is carried out during thesimulation to check if there are any irregularities. Also preliminary tests were conducted to ascertain thatthe collision detection function has no defects.

The first approach taken to verify the algorithm was to create an artificial penetration to see howit behaved under these circumstances. To construct this problem, a single crystal is created and fixed,meaning that it will not be affected by any forces and it will not be considered for the integration ofpositions step. Also, a sphere is created as a local moving body. The sphere is created in such a way thatthere is already an inter-penetration between the sphere and crystal body. This inter-penetration will givea negative penetration distance. This will result in a response that the sphere will be pushed out of thecrystal body.

As seen from the multiple test scenarios shown here, the momentum is shown to act perpendicular tothe contact point or the contact plane as shown in the figures below.

Figure 5.1: The sphere is generated such that it is in contact with the crystal body on one of the faces at t=0with very little inter-penetration.

31

Figure 5.2: In the next time-step, the impulse acts on the sphere in such a way that it moves away from thecrystal body. THe sphere is seen moving away from the crystal.

Figure 5.3: This depicts the parameters that describes the contact. Contact plane and contact normal areperpendicular to each other.

Figure 5.4: The sphere is generated such that it interpenetrates the crystal body at the vertex at t=0, suchthat it overlaps with multiple triangles that describes the crystal body.

32

Figure 5.5: In the next timestep, the sphere is moves in the direction of the contact normal.

Figure 5.6: This depicts the parameters that describes the contact. Contact plane and contact normal areperpendicular to each other.

The next step is to check a single crystal with multiple spheres and making sure that all the spheresundergo proper collision under the given conditions. For this case, again a single crystal is created, andspheres are dropped from the top. All the spheres are driven by the gravitational force, and the results canbe seen from the Figures shown below.

33

Figure 5.7: (a) Multiple spheres are created and dropped from the top. Only gravitational force is actingon the spheres. The crystal body is fixed in the middle of the domain. The ratio of the thickness of thecrystal to the radius of the sphere is 0.3. (b) This figure shows how the spheres that come in contact withthe crystals are bouncing off the face of the crystals. (c) After all the spheres have passed the crystal, onlyspheres that were in the path of the crystal were deflected.

More spheres are dropped on multiple crystals and it is visually checked that there are no spheres thatpass through the crystal. This is shown in the following

Figure 5.8: Here the multiple spheres are created and dropped from the top. Only gravitational force isacting on the spheres. The crystal bodies are fixed in the middle of the domain.

34

Figure 5.9: This figure shows how the spheres that come in contact with the crystals are bouncing off theface of the crystals.

Figure 5.10: After all the spheres have passed the crystal except the spheres that were in the path of thecrystal

5.2 Results with pe-waLBerla CouplingOnce the simulation involving particles and filter crystals was validated, they are imported in the

pe-waLBerla coupling. The crystal filters are mapped to the lattice domain. The flag field looks like asshown in the Figure 5.11

Figure 5.11: Flag field representing the crystals mapped to the grid. Red cells mark the crystals and bluemark the fluid cells.

Initially, fluid simulation is carried out only with the crystal body type. The domain is fixed at 200×400×200 cells. The filter crystals are scaled with the suitable scaling parameter and moved to the center.

35

Then the boundary conditions are set with inflow velocity and pressure outflow with no-slip walls. Theflags of cells which are mapped as crystals are set to no slip as well. After 50000 time-steps, with eachtime-step equalling 2.235e−07s of physical time, the results are shown below. The velocity plot is shown inthe Figure 5.12.

Figure 5.12: Velocity profile after 50000 timesteps. Velocity increases when it passes the filter because itgets accelerated.

Particles are dispersed in the fluid domain and the fluid is driven at a constant velocity. The filtersare located in the center and help in trapping particles. To study the effect of the filter, two indicatorsare crucial. They are the filtration efficiency and the pressure drop. A good filtration performance can becharacterized by a low-pressure drop and high filtration efficiency Li et al. (2016).

Before looking at all the simulations, some of the important parameters are explained here. Thesimulations involve, setting the scaling parameter so that the spatial discretization ∆x is 1e−6m. Thismeans the size of the domain will be as shown in the Figure 4.1. Then it is followed by setting the radiusof particles. At the beginning, radius of 5e−6m was used. Particles are huge compared to the pore sizes,which affects the filtration efficiency. Another parameter that is controlled to depict the dependencies ofthe filtration efficiency is the solid volume fraction (SVF). It can be defined as

Sol i d V ol ume F r acti on = Number o f F i ber Cel l s

Tot al number o f Cel l s(5.1)

where the cells are counted within the width of the filter. This SVF will control the porosity of the filter andhow easily particles bypass the filter. The filtration efficiency can be defined as

F i l tr ati on E f f i ci enc y = Number o f Tr apped par ti cles

Tot al number o f par ti cl es(5.2)

Filtration efficiency is measured for different setups and the results are plotted. The domain andparticles size remain constant but the number of crystals used was varied.

36

Filtration efficiency at various SVF

Filt

rati

on

effi

cien

cy(%

)

20

30

40

50

60

70

80

90

100

110

Solid volume fraction0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08

Figure 5.13: Filtration Efficiency vs SVF for particle sized 10µm

The graph above shows the curve between the filtration efficiency and SVF as defined above. Particlesused were of diameter 10µm. As the SVF increase more particles were trapped. For the SVF value 0.078,efficienncy is 100%. Figure 5.14 and 5.15 shows the two setups where the maximum and the minimumefficiency was observed respectively.

Figure 5.14: Particles being trapped in the filter, SVF = 0.0787, dpar t = 10µm, and Efficiency = 100%.

Figure 5.15: Particles being trapped in the filter, SVF = 0.0183, dpar t = 10µm, and Efficiency = 22.22%.

This next graph shows the efficiency vs. the SVF for particles with diameter of 8µm. Despite with thesame SVF, more particles escape the filter. The maximum efficiency for the same SVF value of 0.078 is

37

97.78%.

Filtration efficiency at various SVF

Filt

rati

on

effi

cien

cy(%

)

20

30

40

50

60

70

80

90

100

Solid Volume Fraction0.02 0.03 0.04 0.05 0.06 0.07 0.08

Figure 5.16: Filtration Efficiency vs SVF for particles sized 8µm.

The next step is to see what happens if the diameter of particles is changed, and the ratio betweenthe diameter of particles and the thickness of the fiber is changed. The crystals are scaled using differentscaling factors and the difference in size is shown in the figure 5.17. The filter created with smallest ∆xis in the front and it increases in size with incremental ∆x. Modifying the discretization, improves thenumber of crystals getting mapped to the domain therefore increasing the SVF.

Figure 5.17: Comparison of Filter sizes for different scaling factors. White = ∆x = 1e−6m , Red = ∆x =7.5e−7m, and Blue = ∆x = 5e−7m.

The table shows the different discretization resolution that was used. The scaling factors change basedon the resolution along with the particle diameters. The number of cells describing particles are fixed atfour cells and five cells.

The change in SVF based on the discretization size is shown in the Figure 5.18. So with the given ∆xvalues, the filtration efficiency is checked.

38

Table 5.1: Simulation setup with different resolution

∆x (10−6m) Domain Size (Cells) Particle Diameter (10−6m)

Setup -1 1 200×400×200 10Setup -2 1 200×400×200 8Setup -3 0.75 300×600×300 7.5Setup -4 0.75 300×600×300 6Setup -5 0.5 400×800×400 5Setup -6 0.5 400×800×400 4

Solid Volume Fraction based on discretization

Soli

dVo

lum

eFr

acti

on

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.1

∆x (m)4e−07 5e−07 6e−07 7e−07 8e−07 9e−07 1e−06 1.1e−06

Figure 5.18: Solid Volume Fraction Change with different ∆x

The efficiency depends on the size of particles and also relative to the SVF. It increases with change inthe size of the domain. This increase improves the filtration efficiency even though particles becomessmaller than the corresponding previous simulation.

Filtration efficiency with different diameters

Filt

rati

on

Effi

cien

cy(%

)

20

30

40

50

60

70

80

Solid Volume Fraction0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1

8 cells10 cells

Figure 5.19: Filtraton Efficiency based on different discretization size of particles and domain

Increase in porosity changes the number of particles that get trapped in the filter. However, in a typicalfilter, there are multiple layers so that the particle gets trapped in either of those layers. With the addition

39

of each layer the domain size is increased in the length direction. It can shown in Figures 5.20 , 5.21 and5.22. The crystals are created back to back but in a random fashion thereby improving the capturingcapability withour increasing the SVF.

Figure 5.20: Domain with one layer of filter, SVF = 0.01828, and domain size = 200×400×200 cells.

Figure 5.21: Domain with two layers of filter, SVF = 0.02848, and domain size = 200×600×200 cells.

Figure 5.22: Domain with three layers of filter, SVF = 0.02707, and domain size = 200×800×200 cells.

The simulations are carried out with particles size of 8µm and 10µm. The Efficiency improves as morelayers are added. This is shown in the Figure 5.23.

40

Filtration efficiency with different layers

Filt

rati

on

Effi

cien

cy(%

)

20

30

40

50

60

70

80

90

Number of Layers0.5 1 1.5 2 2.5 3 3.5

8mum10mum

Figure 5.23: Filtration Efficiency for multiple layers of crystals with particles dpar t = 8µm and dpar t =10µm

Snapshots of particles being captured in the filter crystals are shown in Figures 5.24, 5.25 and 5.26.Particles are captured by all the layers. This explains the increase in efficiency with the multiple layers.

Figure 5.24: Domain with one layer of filter, SVF = 0.018, domain size = 200×800×200 cells , dpar t = 10µmand Filter Efficiency = 22.2%.

41

Figure 5.25: Domain with two layers of filter, SVF = 0.028, domain size = 200×800×200 cells, dpar t = 10µmand Filter Efficiency = 50.7%.

Figure 5.26: Domain with three layers of filter, SVF = 0.027, domain size = 200×800×200 cells,dpar t = 10µmand Filter Efficiency = 65.3%.

42

6 Conclusion

6.1 SummaryIn this work, a new complex geometry, known as the “CrystalBody” was implemented in pe. This rigid

body was described as a convex polyhedron with triangles on the surface. The design of this rigid bodywas based on the filter structures that are used in the acicular mullite diesel particle filter (DPF).

A basic overview of all the concepts was presented. Also presented were the software frameworkpe and waLBerla. The coupling between these two to solve the particulate flows was discussed. Theimplementation of the crystals as rigid bodies in pe was proposed. The description of the filter crystalsusing triangles and how it was imported in the framework was introduced. Then the class hierarchies usedin the parsing of the STL files that describes the crystals were reviewed.

The next step was to construct and implement the collision detection routine between the newlydesigned rigid body and a sphere. Then mapping of the crystals was explained. This process was carriedout using the algorithm presented which ordinarily deals with a variation of ray-triangle intersection thatis regularly utilized in the game engine systems. Subsequently, the sizing of the rigid bodies was presented.The requirement to scale the crystal body and the impact on the domain size was discussed along with theboundary conditions used in the simulation.

Following the setup, simulations were carried out showing the various tests including the test to checkfor the proper response. It validated the contact creation and response routine that was designed. Thenmore particles were involved in the testing, and the response was recorded. Continuing with the testing,the filter crystals are mapped to the LBM grid, and particles are dispersed. Filtration efficiency was studiedby varying the Solid Volume Fraction (SVF). It showed that increasing the number of crystals influencesthe number of particles getting trapped. Then the influence on SVF by scaling factor was investigated. Bychanging the resolution of discretization, efficiency was improved and more crystals were mapped. Nextthe number of layers of the filter wall was incremented and the surge in the efficiency was studied.

6.2 OutlookThe exact representation of the DPF was not possible, but the implementation used to create the filter

wall was a fair approximation. However, the aim of the thesis was to implement a convex body in pe thatwas constructed using triangles, so that more and more rigid bodies with shapes that cannot be describedby simple geometric primitives can be utilized in the simulation. The next step can be to design a systemso that any convex rigid body using triangles would be directly imported. Collision between convex bodypairs can implemented on top of the convex body-sphere collision module that is used in this work. Thiscan be useful in improving the functionality of pe.

In another direction, this can be developed further by incorporating the exact filter design in thesimulation and maintain particles relative to the actual length sizes that are seen in the diesel particulatesystems. The current implementation has restriction on domain sizes because resolution of spatialdiscretization has a huge impact on the computation costs. This restriction should be overcome by using asmaller domain or taking a representative part of the filter to study the effect on filtration efficiency. Thenthe effect of the parameters such as the porosity, the permeability and the fiber thickness on the filtrationefficiency and the pressure drop can also be studied.

43

7 ReferencesMichael Fiebig, Andreas Wiartalla, Bastian Holderbaum, and Sebastian Kiesow. Particulate emissions

from diesel engines: correlation between engine technology and emissions. Journal of OccupationalMedicine and Toxicology, 9(1):6, Mar 2014. ISSN 1745-6673. doi: 10.1186/1745-6673-9-6. URL https://doi.org/10.1186/1745-6673-9-6.

Tao Chen, Zhixin Wu, Jinke Gong, and Jia-qiang E. Numerical simulation of diesel particulate filterregeneration considering ash deposit. Flow, Turbulence and Combustion, 97(3):849–864, Oct 2016. ISSN1573-1987. doi: 10.1007/s10494-016-9717-6. URL https://doi.org/10.1007/s10494-016-9717-6.

Samir Bensaid, Daniele L. Marchisio, and Debora Fino. Numerical simulation of soot filtration andcombustion within diesel particulate filters. Chemical Engineering Science, 65(1):357 – 363, 2010. ISSN0009-2509. doi: https://doi.org/10.1016/j.ces.2009.06.051. URL http://www.sciencedirect.com/science/article/pii/S0009250909004473. 20th International Symposium in Chemical ReactionEngineering—Green Chemical Reaction Engineering for a Sustainable Future.

Chao Zhu, Chao-Hsin Lin, and Chun Shun Cheung. Inertial impaction-dominated fibrous filtration withrectangular or cylindrical fibers. Powder Technology, 112(1):149 – 162, 2000. ISSN 0032-5910. doi:https://doi.org/10.1016/S0032-5910(99)00315-0. URL http://www.sciencedirect.com/science/article/pii/S0032591099003150.

Antonis Karadimos and Raffaella Ocone. The effect of the flow field recalculation on fibrous filter loading: anumerical simulation. Powder Technology, 137(3):109 – 119, 2003. ISSN 0032-5910. doi: https://doi.org/10.1016/S0032-5910(03)00132-3. URL http://www.sciencedirect.com/science/article/pii/S0032591003001323.

Rafał Przekop, Krzysztof Grzybowski, and Leon Gradon. Energy-balanced oscillatory model for de-scription of particles deposition and re-entrainment on fiber collector. Aerosol Science and Technol-ogy, 38(4):330–337, 2004. doi: 10.1080/02786820490427669. URL http://dx.doi.org/10.1080/02786820490427669.

S.J. Dunnett and C.F. Clement. A numerical study of the effects of loading from diffusive depositionon the efficiency of fibrous filters. Journal of Aerosol Science, 37(9):1116 – 1139, 2006. ISSN 0021-8502. doi: https://doi.org/10.1016/j.jaerosci.2005.08.001. URL http://www.sciencedirect.com/science/article/pii/S0021850205001679.

Thilo Müller, Jörg Meyer, and Gerhard Kasper. Low reynolds number drag and particle collision efficiency ofa cylindrical fiber within a parallel array. Journal of Aerosol Science, 77(Supplement C):50 – 66, 2014. ISSN0021-8502. doi: https://doi.org/10.1016/j.jaerosci.2014.07.007. URL http://www.sciencedirect.com/science/article/pii/S0021850214001153.

S.A. Hosseini and H. Vahedi Tafreshi. Modeling particle-loaded single fiber efficiency and fiber drag usingansys–fluent cfd code. Computers & Fluids, 66(Supplement C):157 – 166, 2012. ISSN 0045-7930. doi:https://doi.org/10.1016/j.compfluid.2012.06.017. URL http://www.sciencedirect.com/science/article/pii/S004579301200240X.

Fuping Qian, Naijin Huang, Xiaojie Zhu, and Jinli Lu. Numerical study of the gas–solid flow characteristic offibrous media based on sem using cfd–dem. Powder Technology, 249(Supplement C):63 – 70, 2013. ISSN0032-5910. doi: https://doi.org/10.1016/j.powtec.2013.07.030. URL http://www.sciencedirect.com/science/article/pii/S0032591013004890.

Liu-Chao Qiu. A coupling model of dem and lbm for fluid flow through porous media. Procedia Engi-neering, 102(Supplement C):1520 – 1525, 2015. ISSN 1877-7058. doi: https://doi.org/10.1016/j.proeng.2015.01.286. URL http://www.sciencedirect.com/science/article/pii/S1877705815003057.New Paradigm of Particle Science and Technology Proceedings of The 7th World Congress on ParticleTechnology.

Klaus Iglberger. Software design of a massively parallel rigid body framework. 2010.

Mario Heene. Extension of the pe physics engine by discrete element methods. Bachelor thesis, Lehrstuhlfür Informatik, 10, 2011.

44

Tobias Preclik. Models and algorithms for ultrascale simulations of non-smooth granular dynamics. 2014.

Christoph Rettinger and Ulrich Rüde. A comparative study of fluid-particle coupling methods for fullyresolved lattice boltzmann simulations. Computers & Fluids, 2017.

David Baraff. Physically based modeling: Rigid body simulation. SIGGRAPH Course Notes, ACM SIGGRAPH,2(1):2–1, 2001.

Ladislav Kavan. Rigid body collision response. Vectors, 1000:2, 2003.

David H Eberly. Game physics. CRC Press, 2010.

Murilo G Coutinho. Guide to dynamic simulations of rigid bodies and particle systems. Springer Science &Business Media, 2012.

Matthew Moore and Jane Wilhelms. Collision detection and response for computer animation. SIGGRAPHComput. Graph., 22(4):289–298, June 1988. ISSN 0097-8930. doi: 10.1145/378456.378528. URL http://doi.acm.org/10.1145/378456.378528.

Stéphane Redon, Abderrahmane Kheddar, and Sabine Coquillart. Fast continuous collision detectionbetween rigid bodies. Computer Graphics Forum, 21(3):279–287, 2002. ISSN 1467-8659. doi: 10.1111/1467-8659.t01-1-00587. URL http://dx.doi.org/10.1111/1467-8659.t01-1-00587.

Gino Johannes Apolonia van den Bergen. Collision detection in interactive 3D computer animation.University Press Facilities, 1999.

Jason Gregory. Game engine architecture. crc Press, 2009.

Christer Ericson. Real-time collision detection. CRC Press, 2004.

Jeffrey C Trinkle, J-S Pang, Sandra Sudarsky, and Grace Lo. On dynamic multi-rigid-body contact problemswith coulomb friction. ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für AngewandteMathematik und Mechanik, 77(4):267–279, 1997.

Bruce R Donald and Dinesh K Pai. On the motion of compliantly-connected rigid bodies in contact. ii. asystem for analyzing designs for assembly. In Robotics and Automation, 1990. Proceedings., 1990 IEEEInternational Conference on, pages 1756–1762. IEEE, 1990.

Peng Song, Peter Kraus, Vijay Kumar, and Pierre Dupont. Analysis of rigid-body dynamic models forsimulation of systems with frictional contacts. TRANSACTIONS-AMERICAN SOCIETY OF MECHANICALENGINEERS JOURNAL OF APPLIED MECHANICS, 68(1):118–128, 2001.

David E Stewart and Jeffrey C Trinkle. An implicit time-stepping scheme for rigid body dynamics withinelastic collisions and coulomb friction. International Journal for Numerical Methods in Engineering,39(15):2673–2691, 1996.

David E Stewart. Rigid-body dynamics with friction and impact. SIAM review, 42(1):3–39, 2000.

Klaus Iglberger and Ulrich Rüde. Massively parallel rigid body dynamics simulations. Computer Science -Research and Development, 23(3):159–167, Jun 2009. ISSN 1865-2042. doi: 10.1007/s00450-009-0066-8.URL https://doi.org/10.1007/s00450-009-0066-8.

Hannes Wengenroth. Rigid body collisions. Computer Science Department, 10, 2007.

Jan Götz. Massively Parallel Direct Numerical Simulation of Particulate Flows. Verlag Dr. Hut, 2012.

Elmer G Gilbert, Daniel W Johnson, and S Sathiya Keerthi. A fast procedure for computing the distancebetween complex objects in three-dimensional space. IEEE Journal on Robotics and Automation, 4(2):193–203, 1988.

Florian Schornbaum. A Real-time Capable Impulse-based Collision Response Algorithm for Rigid BodyDynamics. PhD thesis, Master’s thesis, Friedrich–Alexander–Universität Erlangen-Nürnberg, 2010.

Timm Krüger, Halim Kusumaatmaja, Alexandr Kuzmin, Orest Shardt, Goncalo Silva, and Erlend MagnusViggen. The Lattice Boltzmann Method. Springer, 2017.

45

P. L. Bhatnagar, E. P. Gross, and M. Krook. A model for collision processes in gases. i. small amplitudeprocesses in charged and neutral one-component systems. Phys. Rev., 94:511–525, May 1954. doi:10.1103/PhysRev.94.511. URL https://link.aps.org/doi/10.1103/PhysRev.94.511.

Y. H. Qian, D. D’Humières, and P. Lallemand. Lattice bgk models for navier-stokes equation. EPL (Euro-physics Letters), 17(6):479, 1992. URL http://stacks.iop.org/0295-5075/17/i=6/a=001.

Renwei Mei, Wei Shyy, Dazhi Yu, and Li-Shi Luo. Lattice boltzmann method for 3-d flows with curvedboundary. Journal of Computational Physics, 161(2):680–699, 2000.

Christian Feichtinger. Design and performance evaluation of a software framework for multi-physicssimulations on heterogeneous supercomputers. 2012.

Klaus Iglberger, Nils Thürey, and Ulrich Rüde. Simulation of moving particles in 3d with the latticeboltzmann method. Computers & Mathematics with Applications, 55(7):1461 – 1468, 2008. ISSN 0898-1221. doi: https://doi.org/10.1016/j.camwa.2007.08.022. URL http://www.sciencedirect.com/science/article/pii/S0898122107006256. Mesoscopic Methods in Engineering and Science.

Sauro Succi. The lattice Boltzmann equation: for fluid dynamics and beyond. Oxford university press, 2001.

Irina Ginzburg, Dominique d’Humières, and Alexander Kuzmin. Optimal stability of advection-diffusionlattice boltzmann models with two relaxation times for positive/negative equilibrium. Journal ofStatistical Physics, 139(6):1090–1143, Jun 2010. ISSN 1572-9613. doi: 10.1007/s10955-010-9969-9. URLhttps://doi.org/10.1007/s10955-010-9969-9.

C. Feichtinger, S. Donath, H. Köstler, J. Götz, and U. Rüde. Walberla: Hpc software design for computationalengineering simulations. Journal of Computational Science, 2(2):105 – 112, 2011. ISSN 1877-7503.doi: https://doi.org/10.1016/j.jocs.2011.01.004. URL http://www.sciencedirect.com/science/article/pii/S1877750311000111. Simulation Software for Supercomputers.

Anthony JC Ladd. Numerical simulations of particulate suspensions via a discretized boltzmann equation.part 1. theoretical foundation. Journal of fluid mechanics, 271:285–309, 1994.

Kristina Pickl, Jan Götz, Klaus Iglberger, Jayant Pande, Klaus Mecke, Ana-Suncana Smith, and Ulrich Rüde.All good things come in threes - three beads learn to swim with lattice boltzmann and a rigid body solver.CoRR, abs/1108.0786, 2011. URL http://arxiv.org/abs/1108.0786.

Mario Botsch, Mark Pauly, Christian Rossl, Stephan Bischoff, and Leif Kobbelt. Geometric modeling basedon triangle meshes. In ACM SIGGRAPH 2006 Courses, page 1. ACM, 2006.

Aleksander J Pyzik and Cheng G Li. New design of a ceramic filter for diesel emission control applications.International Journal of Applied Ceramic Technology, 2(6):440–451, 2005.

A Gil, JPG Galache, C Godenschwager, and U Rüde. Optimum configuration for accurate simulations ofchaotic porous media with lattice boltzmann methods considering boundary conditions, lattice spacingand domain size. Computers & Mathematics with Applications, 73(12):2515–2528, 2017.

R Nürnberg. Calculating the volume and centroid of a polyhedron in 3d, 2013.

Florian Schornbaum. Hierarchical hash grids for coarse collision detection. Student Thesis, University ofErlangen-Nuremberg, 2009.

Tobias Scharpff. Numerische simulation scharfkantiger teilchen in der pe physikengine. Diploma Thesis,University of Erlangen-Nuremberg, 2014.

HPC Cluster-LSS information. https://www10.informatik.uni-erlangen.de/de/forschung/hpc-cluster/. Accessed: 2017-10-28.

Wei Li, Shengnan Shen, and Hui Li. Study and optimization of the filtration performance of multi–fiberfilter. Advanced Powder Technology, 27(2):638 – 645, 2016. ISSN 0921-8831. doi: https://doi.org/10.1016/j.apt.2016.02.018. URL http://www.sciencedirect.com/science/article/pii/S0921883116000662.

46