Simulation of filtration processes in deformable media Part 3 . 1 : Basic concepts and particle-fluid force implementation of a non-spherical dirt particle solver

A Lagrangian solver to realistically model large, non-spherical dirt particles and their behaviour in the vicinity of deformable filtration fibres has been programmed. While this paper focuses on basic solver concepts as well as drag force implementations, a related article, concerning the realisation of interaction effects and result verification, is forthcoming, [3]. Within the framework of a digitally reconstructed, deformable filter fibre geometry, the solver traces the governing multi physics effects down to the occurrence of single forceand torque vectors. In order to go from an initial, spherical particle model [2], to a more sophisticated, non-spherical model, the capabilities of a Six Degrees of Freedom Solver have been included in the programming. A panel model and the concept of satellite help points are used to handle particles that encompass several fluid calculation cells. An innovative drag force implementation allows the consideration of rotationaland shear flow effects on particle motion. Results are evaluated and compared to an analytical formulation.


INTRODUCTION
The Open Source Computational Fluid Dynamics (CFD) toolbox OpenFOAM® has served as a programming environment for the development of a novel, deterministic, micro scale, fluid-particle-fibre filtration solver [1], [2].This C++ based simulator is supposed to become an important tool for the CFD design of filter media for automotive lubrication.It was created to consider all physically relevant effects that go along with, or lead to a micro scale, dirt particle deposition in a realistically reconstructed filtration fibre geometry.
Concerning this subject, two papers, have been previously published.While [1] focuses on fibre deformation and fluid-structure interaction (FSI) effects, [2] describes the development of a large, spherical particle model for filtration applications.
A second article [3] on the non-spherical particle model is forthcoming.It is closely related to this work and is to be seen as a sequel.It mainly concerns itself with the handling of particle interaction with its surroundings.Moreover it deals with the creation of a simplified, semi-analytical approach to verify solver functionality and result quality.
The current work presents a significant extension of the spherical dirt particle model formulated in [2].It describes the basic concepts as well as the essential drag force implementation method behind our novel, realistic, Lagrangian, non-spherical particle model.

Simulation of filtration processes in deformable media
Some of the main aspects of the development effort behind this paper are: • Implementation of a Six Degrees of Freedom (DOF) solver for the Lagrangian particle momentum equations (PME) in OpenFOAM®.

•
Design of an explicit, force-and torque vector model, that reduces the modelling to the mere formulation of single force effects.

•
Introduction of an adaptive time stepping scheme for explicit Euler discretization of the PME, [4].

•
Device of a surface help point scheme to account for large particle effects in terms of fluid-particle, particle-fibre and particle-particle interaction.

•
Creation and verification of a drag force implementation that uses a combination of non-spherical, semi-empirical drag force formulas [23], a panel method to consider free flow swirling effects, a plugging method to include inter particle and particle-fibre hydrodynamics and a simple adoption of basic concepts known from the immersed boundary method [16].

FILTRATION SOLVER
The original filtration solver, presented in [1] and [2] is capable of solving complex fluidparticle-fibre interaction problems for microscopic, digitally reconstructed fibre samples ranging up to a few hundred microns in diameter and thickness.To reconstruct the microscopic geometry as accurately as possible, input data from computer tomographic scans is used.Rudimentary stacks of grey scale images are being digitalized by Matlab® utilities.An exemplary result is shown in Figure 1.
Both the fibres and their surrounding space are being transformed into structured grid meshes that represent the solid and fluid framework for Computational Fluid Dynamics (CFD) analysis.Because of the geometry range of ~200µm, the highly viscous, Newtonian oil current featuring a kinematic fluid viscosity of η f ~10 -4 m 2 /s, and relatively slow flow velocities of <0.4 m/s, the local Reynolds numbers in the fibre vicinity are expected to be mostly <0.5, but surely <1.0.Furthermore Knudsen numbers range well below 0.015.Thus continuum equations are valid and the consideration of diffusion effects on particle motion becomes unnecessary.After all a simple, incompressible, laminar and isothermal fluid solver can handle the situation.
For quite some time, deformation effects have been suspected to have significant impact on the filter characteristics of a fibre.Therefore the original solver was designed to handle fluid-structure interaction effects (FSI).The FSI utility features a stiff, explicit coupling and uses only one, single finite volume solver to handle the governing fluid dynamics as well as the structural mechanics on the solid side, [1].
The original particle model was a spherical, Lagrangian, fully deterministic (nonstochastic) approach with the capability to interact with the surrounding, Eulerian fluid-fibre framework.Each particle can extend well beyond the borders of a single calculation cell and can sense and affect fluid conditions within a multiple cell region of the fluid mesh.In extension of the spherical particle model, a highly detailed, more sophisticated and more accurate, non-spherical particle model was developed and can hereby be presented.

DRAG FORCES AND PARTICLE RELAXATION TIMES
An obvious reason to go from a mere spherical dirt particle description to a more realistic, non-spherical approach lies within a significant difference in drag-force-to-mass-ratio.A good way to demonstrate the difference is to take a look at spherical and non-spherical particle relaxation times τ p of mass equivalent particles.With m p being the particle mass and ρ p being the particle density, the diameter of a mass equivalent sphere D sph can be written as: (1) Since the particle Reynolds numbers under consideration, range significantly below 1, Stokes drag conditions can be assumed.Thus the expression for the particle relaxation time for spherical particles τ p,sph in the flow domain is given by: (2) Here µ f is the dynamic fluid viscosity.
For non-spherical particles a drag force correlation, proposed by Hölzer & Sommerfeld, [23] shall be chosen.In this case the definition of the non-spherical particle relaxation time τ p,nonsph is much more complex and reads: Here u rel is the relative fluid-particle velocity and the constants C 0 , C 1 , C 2 and C 3 are: (4) Where ρ f is the fluid density and A f,ell is the frontal area of an ellipsoid particle, projected onto a plane, perpendicular to the relative fluid-particle velocity vector.
The comparison of Equ.2 and Equ.3 shows that non-spherical particle relaxation times are generally lower than those of mass equivalent spheres.In [4] the parameter α ax is introduced to measure deviation from spherical shape.It represents the medium, relative half axis deviation around D sph and is defined as: (8) Here a, b and c are the lengths of the three particle half axis, whereby a ≥ b ≥ c.Using α ax as a parameter, it becomes apparent that, the further the particle shape deviates from being a sphere (higher α ax ), the smaller τ p,nonsph will be, [4].A corresponding plot of the situation, as seen in Figure 2, reveals that non-spherical particle relaxation times show a dependency on local fluid conditions, while spherical particle relaxation times do not.
Furthermore the results in Figure 2 show that non-spherical particle relaxation times for highly non-spherical particles (α ax ≥ 1) amount to less than 1/5 th of spherical relaxation times.Therefore a mere spherical particle model would significantly underestimate fluid skin friction-and form drag forces on supposedly arbitrarily shaped dirt particles.One obvious consequence of disregarding particle shape effects for filtration simulation would be a certain overestimation of filter fibre efficiencies

MODEL BASICS
This chapter lists some basic concepts and innovative implementation schemes that had to be chosen and/or developed in order to create a suitable framework for the non-spherical particle model.

PARTICLE GEOMETRY: ELLIPSOID SHAPE
The non-spherical particle shape representation is chosen to be an ellipsoid with three independent, geometrical degrees of freedom.This choice offers the versatility to approximate many shapes from sticks to plates (see Figure 3) on the one hand and can be mathematically described pretty easily, as seen in Equ.9, on the other hand.(9) Here x', y' and z' are coordinates of a Lagrangian, co-rotational coordinate system, with base vectors being aligned along the particle's main axis (see chapter 3.2).
x a  Figure 3 The ellipsoid shape can approximate a wide variety of geometries, e.g.plates and sticks.
The ellipsoid particle volume V p is given by: (10) An essential quantity for calculating skin friction forces on the particle, is the particle surface area A p .In this work A p for ellipsoids shall be approximated by a comparatively simple approximation formula, proposed by Thomsen, [5]: (11) With p ≈ 1.6075, this formula is reported to yield a maximum of +/-1.061% deviation about the correct result.

EULER AND LAGRANGE COORDINATE SYSTEM
The fluid and FSI calculations are based upon well known, Eulerian principles and require only one Cartesian coordinate system, with base vectors e x , e y and e z and coordinates, x, y and z.In the course of the FSI calculation however, the fluid mesh actually works as Lagrangian mesh that adjusts itself to displaced fibres.Whereas the fibre structure mesh itself remains in its original position, after the fashion of a typical Eulerian mesh, [1].
For the particle calculation the partly Lagrangian character of the fluid mesh is completely irrelevant.The particle solver does not require separate meshing.To account for particle position X p and orientation, an additional co-rotational coordinate system is introduced.The particle coordinate system, with base vectors e px , e py and e pz being aligned along the main particle axis, as seen in Figure 4, originates from the particle mass centre.Its coordinates are written as x', y' and z'.The relationship of any single point P within the Eulerian system, to the corresponding point P' within the co-rotational Lagrangian system, is given by the following formula: Here the index n denotes the axe directions x, y and z respectively and the base vectors of the particle coordinate system are given by: (13) According to Equ.12 the transformation operation T can be defined as: (14) With A being the transformation matrix: (15) Accordingly the retransformation from P' to P is computed as: (16) Here the retransformation operation T' is formulated by using the transposed transformation matrix A T : (17) A similar multiple coordinate system approach is used by Rosdahl, [8].His solver uses a third, additional, co-moving coordinate system, which also originates from the mass centre of the particle and is aligned along the basis of the outer, inertial, Eulerian coordinate system.

SIX DEGREES OF FREEDOM SOLVER
Non-spherical particle motion in general consists of three degrees of translational as well as three degrees of rotational freedom.This is why the original, spherical particle solver had to be transformed into a more general, six DOF solver.It can now solve the translational and rotational, Lagrangian equations of motion for N external forces F j , which act on the particle: (18) Here u p is the particle velocity and ρ p the particle density.(19) Hereby I n,j is the particle moment of inertia around the rotational axe n j of any acting torque vector.The direction of the lever is given by r j .I n,j is calculated by use of the particle 1 2 e e z ',3 moment of inertia tensor I and by use of Equ.21.I in its most general form can be described by using a continuous, spatial density ρ sp (x, y, z), [8]: (20) Here V sp describes the local space which completely encompasses the object, r k is the distance vector to the axe of rotation, E stands for the identity matrix and denotes the outer product.(21) Inserting the expressions for standard, ellipsoid, principle moments of inertia, this finally amounts to:

SINGLE FORCE AND TORQUE VECTORS
A number of individual, translational and rotational force-and torque effects with influence on particle trajectory and deposition behaviour have to be accounted for.Those effects can be parted into three basic categories: particle-fluid (see chapter 4), particle-fibre [3] and particle-particle [3] interactions.
It is inaccurate to traditionally model one individual particle-deposition effect, without taking into account the interaction with other particles, or a changing fluid field [11], [12].Thus in this work all modelling is broken down to the level of individual force effects and their resulting torques.The following interaction forces are relevant: • particle-wall impact force, F wall • particle-fibre interaction force, F fibre • particle-particle impact force, F collision • particle-fluid interaction (drag) force, F fluid • force due to pressure gradient (form drag), F pressure • force due to shear flow (shear drag), F shear • gravity, F g A simultaneous calculation of F fluid , F pressure and F shear would yield an overestimation of fluid-particle interaction forces.Depending on the specific mode of operation, either F fluid or F pressure /F shear are calculated.
In order to model resulting torques T j , the exact positions r j of acting forces F j need to be known.This is why a surface help point method (see chapter 2.5) was devised.Figure 5 illustrates an assembly of small, non-spherical particles and the corresponding system of acting forces and torques, that cause translation and rotation.

PARTICLE SHAPE CONCEPTS
In order to consider rotational effects, collision-impact scenarios or any shape-related phenomenon, the moving object has to extend beyond a simple, point-like representation.Thus the surface help point method, as well as a simple panel method to discretize the particle surface are introduced.

Surface and pressure/velocity help points
A cloud of up to M= 68 help points per particle is used (see Figure 6).18 surface help points are positioned directly at the surface of the particle to serve as collision detectors and pressure/velocity probes.An additional 48 pressure/velocity help points are located at crucial positions of the particle-panel model (see chapter 2.5.2) and detect local fluid field conditions.The help points surround the ellipsoid at constant positions HP m ' within the framework of the co-rotational particle coordinate system.Thus each help point conserves its relative position to the particle centre, while the particle moves arbitrarily within the Eulerian fluid domain, "dragging along" the Eulerian help point positions HP m .Figure 6 Non-spherical particle with 18 surface help points and 48 pressure/ velocity help points.
The help points essentially serve two purposes: 1.
In their function as pressure/velocity help points they detect local fluid conditions.2.
They predefine the current particle movement by tracking the individual, projected trajectory.Collision that might occur along this course are detected.
Using a simple, temporal Euler discretization [4], the help point position HP m i at time i can be projected to its new position HP m i+1 , at time i+1, after particle time step ∆t p .The new position can be calculated as: (23) Here r hp is the help point distance vector to the particle mass centre.This particle progression scheme is only used if collision events are to be expected.The linear trajectory is probed for obstacles.If a collision occurs at position X coll before HP m i+1 is reached, the fraction f m is set to: Then the new particle time step ∆t f * is calculated as: (25) Now the actual particle movement is conducted.
If no collision events are to be expected, the solver conducts translational and rotational operations on the particle mass centre X p and on the particle base vectors e p,n .To find the new help point positions at time i+1, a simple coordinate transformation suffices: (26) Note the fact that the co-rotational help point positions HP m ' remain unchanged at all times.

Panel method
While the surface help point scheme has been designed to aid in the modelling of collisions and in the detection of local flow field conditions, a simple panel method is introduced to get a hold of hydrodynamic drag-and lift forces.Within the co-rotational coordinate system, the fixed help point positions are used as a framework to encase the ellipsoid with a system of edges and panels (see Figure 7).

ADAPTIVE TIME STEPPING
Because of the simple, explicit, temporal Euler discretization that is used for particle progression, certain numerical instabilities have to be opposed [18]- [21].Studies on the subject have led to the development of an adaptive time stepping scheme, that compensates individual particle needs for higher time resolution.With σ rel,UD being the user-definable, relative standard deviation between a numerical and an analytical particle speed up curve, the number of particle sub time steps J can be calculated for two levels of accuracy.For σ rel,UD ≤ 1.2% u f : For 1.2% u f < σ rel,UD ≤ 4.0% u f : (28) Hereby u f stands for the uniform velocity of the surrounding fluid field.A full report on the results has been published in [4].

PARTICLE MOMENTUM EQUATION
All changes in translational and rotational motion of a particle stem from the sum of acting forces.Applying Newton's second law, the translational particle momentum equation (PME) for arbitrarily shaped particles and arbitrary flow conditions, is given.The PME presents the framework for all particle-motion modelling behind the non-spherical particle solver.It can be written as: (29)

Steady State Forces
Unsteady Forces Event Forces Here F d is the drag force, F h the hydrodynamic lift force, F Magnus is the Magnus force, F Saffman is the Saffman force, F Faxen is the Faxen force, F g is the gravity force, F b is the buoyancy force, F VM is the added (virtual) mass effect, F Basset is the Basset (history) force and F e,i stands for the i-th of N event forces [6], [14], [15].The individual force contributions, summarized in Equ.29 can be divided into three main categories: steady state forces, unsteady forces and event forces.
The specialization of Equ.29 for small, spherical particles, that are not two-way-coupled to a uniform, surrounding fluid, gives the Basset Bousinesque Oseen (BBO) equation [13], [27] without Faxen terms or interaction with solids or other particles.In this case the individual force contributions can be formulated as seen in Figure 9.
Figure 9 Formulation of force contributions as they would look like in order to turn the PME in Equ.29 into the classic BBO equation for small, spherical, non-coupled particles.
In Figure 9, c d denotes the drag coefficient and Re p stands for the particle Reynolds number.
Usually a PME formulation like the one in Equ.29 is used for small particles and the classical Euler-Lagrange approach [24].This work however, treats large particles that span multiple fluid cells, and still retains the typical Euler-Lagrange methodology.Thus a specifically adjusted, numerical scheme to model particle-fluid interaction becomes necessary.
Comparable programs, like that of Schütz [7], are using particle-related re-meshing of the fluid grid.

PARTICLE-FLUID INTERACTION
To maximize calculation efficiency, a detailed drag implementation, specifically adapted to the case of dirt particle filtration in lubricants has been created.The particle-fluid interaction model, consists of two alternative modules: Free Flow particle-fluid module Fibre vicinity particle-fluid module Dirt particles are injected into the free flow regime upstream of the filter fibre geometry, where they occur in very low volume fractions.Particle-particle interaction and hydrodynamic particle impact on the fluid can be neglected here.As soon as the particles reach the fibre vicinity, the two-way coupling takes effect and inter-particle-as well as full particle-fluid interaction becomes relevant.Those fundamentally different situations require separate drag modelling schemes in order to guarantee a good balance between accuracy and efficiency.

FREE FLOW PARTICLE-FLUID INTERACTION MODULE
Within the free flow regime, all particle interactions with their surroundings are handled by the free flow module.In this zone, the most important aspects of the prevailing hydrodynamic situation are: 1.) The ratio between particle diameter D sph and minimal distance to the nearest fibre (wall) boundary patch h w can be considered as small.Wall proximity has no effect on particle drag.2.) Due to very low particle volume fractions, the ratio between D sph and the minimal distance between neighbouring particles h p can be considered as small.No physical nor hydrodynamic particle interaction takes place.3.) Due to very low particle Reynolds number Re p and very low particle volume fraction, the hydrodynamic particle effect on the fluid can be neglected.No twoway-coupling is necessary.
In the free flow regime it is primarily important to grasp torques, acting on the particle.Rotational effects due to non-uniform flow fields can lead to a pre-alignment of the particles, so that average penetration depth and filter fibre efficiencies are being influenced.In that context the panel description (see chapter 2.5.2) of the ellipsoid shape is of special importance.The particle is enclosed by M panels and each panel j is subject to drag forces F d,j (which consist of pressure-and shear flow contribution, F p,j and F t,j respectively) and hydrodynamic lift forces F h,j .Forces that are better calculated by considering the entire particle are: gravity F g , buoyancy F b and N event forces ΣF e,i .Thus the adapted PME within the free flow regime looks like Equ.30.In analogy, torque effects on non-spherical rotation are described by: (31) Here r j stands for the distance vector of each surface panel centre HP j to the particle mass centre X p and r i denotes the distance vector from X p to any particle help point HP i that senses an impact event.
The Basset history force and the added, virtual mass are being neglected because of the lack of strongly in stationary, relative particle-fluid flow.Comparing Equ.30 with Equ.29, the following parallels can be drawn: Figure 10 shows a sketch of how the individual force contributions act on each panel and affect the particle.In Figure 10, F in stands for the force contribution of the incoming current, while F out is the force contribution of the outgoing current as it would look like if it were deviated by the panel surface.The total force acting on each panel F panel is given by: (33) The following sub chapters describe the procedural calculation of free flow drag-and lift force as well as torque effects within the module.

FORCE CALCULATION
First the drag force contribution F d,j on each panel j has to be calculated.The drag force term consists of a form drag-and a shear drag contribution, F p,j and F T,j .While F T,j acts perpendicular to the panel surface normal n p,j , F p,j acts parallel to n p,j .However, since F d,i is supposed to act in the direction of u rel,j , the following ratio has to hold: Here e urel,j is the base vector of the relative fluid-particle-panel velocity encountered at the panel centre.The total panel drag coefficient c d,panel depends on the form drag coefficient c d,p and on the shear drag coefficient c d,shear and is given by: (35) In Equ.36 to Equ.40, form drag-and shear drag vectors are listed, as well as the total panel drag vector and its dependence on form and shear contribution, F p,j * and F T,j * respectively.Figure 11 shows a sketch of the situation.In Equ.36 through Equ.40,A j is the panel surface area and u rel,j is the relative fluidparticle-panel velocity u f -u p,j .The particle-panel velocity is given by the velocity of the particle mass centre u p and the rotational velocity contribution: The panel Reynolds number is written as Re j and is defined by using the hydraulic diameter d h,j of the panel and the kinematic fluid viscosity η f : (42) Re Secondly the hydrodynamic lift force F h,j , which stems from the deviation of the fluid at the panel, is calculated.The hydrodynamic lift is connected to F in,j , F out,j and F d,j via a simple, local force balance (see Figure 10): (43) Note that the superscript ' denotes the fact that force values are scaled by the acting surface area A j .In addition to Equ.43, F h,j is defined to act perpendicular to F d,i , so that: The drag F d,j ' is given by Equ.40, while F in ' can be easily derived out of the local fluid field information, obtained by the pressure/velocity help points.
(45) While the value of F out ' is not known in advance, its base vector e out,j is given because of panel orientation n p,j and relative panel-fluid velocity u rel,j .
(46) Equ.43 and Equ.44 constitute a system of four equations and four unknowns: The three components of hydrodynamic lift F h,j,x , F h,j,y, F h,j,z and the absolute value of deviated flow momentum |F out,j |.The solution yields the following expressions for the local, hydrodynamic lift force vector F h,j ' and the vector of deviated fluid momentum F out,j ': ( ) Due to non-coupling, the wake of the particle is not simulated in the free flow module.Therefore a panel has to face the current in order to yield acceptable F d,j and F h,j results.The condition for calculating the individual force balance and for considering the panel is: Consequentially the overall, unscaled drag force F d unsc and hydrodynamic lift force F h unsc are given by the contributions of all N considered panels:

Weighing method and torque effect calculation
The F d,j and F h,j calculation serves a useful purpose: to obtain an idea of the force distribution over the particle surface.This is necessary to grasp rotational fluid field effects on the aligning particle.
In order to improve the quantitative estimate on overall drag and lift force contributions, the panel results are being scaled to fit a newly found, empirical drag law for non-spherical particles that has been presented by Hölzer & Sommerfeld,[23]: Here Φ, Φ length and Φ cross , are the standard-, the lengthwise-and crosswise-sphericity, respectively.Using the Hölzer-Sommerfeld approach, the over all scaled drag force F d sc on the particle can be calculated.In order to get a hold of realistic, rotational torque effects, the originally calculated, unscaled force contributions are scaled by the ratio F d sc / F d unsc : (53) (54) Then the more accurate, scaled, rotational torque effects are computed using Equ.31: (55) Figure 12 shows a test case, where a longish, non-spherical particle approaches a valve of higher flow velocities and lower pressure.As physically plausible and expected, the given drag implementation models the occurring shear flow and pressure gradients over the particle surface in such a way that the particle aligns itself along the fluid stream lines.The translational and angular velocity vectors adapt to the local fluid field conditions which leads to a particle-slip effect (see Figure 12).This behaviour causes the non-spherical slip effect with relevance for filtration efficiency and particle penetration depth [3].

Particle in fibre vicinity
As soon as the particle enters into the vicinity of the fibre geometry, the hydrodynamic situation changes as a whole and the fibre vicinity drag module takes over.The situation features the following characteristics: 1.) Particle-wall flow effects can no longer be neglected since the ratio between particle diameter and minimal particle-wall (fibre) distance h p is per definition no longer small.2.) Particles accumulate at the fibre in considerable volume fractions and the ratio between particle diameter and medium, minimal particle-particle distance h pp is no longer small either.Particles interact hydrodynamically and physically by plugging each others flow path.3.) The high particle volume fractions, lead to a plugging of the fluid flow path, diverting the flow and causing increased pressure drop over the filter.Hydrodynamic particle impact on the fluid (two-way coupling) becomes essential.
Empirical expressions, describing every one of those effects individually, can be found in literature [27], [13].
The fibre vicinity drag module incorporates all those effects in a dynamic combination to a highly complex, multi-parameter interaction situation.Local fluid cells, encompassed by particles are plugged by porosity adjustment.The procedure is analogous to the one presented in [2] but has been extended to non-spherical particle shapes, refined and quantified as shall be seen in the following.
The plugging causes the fluid to be diverted around the fluid cell which leads to a local pressure build up p i , that can be sensed by any of the N pressure help points HP i at the particle surface (see Figure 13).Since each pressure help point represents 1/N th of the entire particle surface area A p and since pressure always acts perpendicular to the local surface normal n i , the total pressure force F p on the particle can be written as: (56) For infinitesimally fine grid spacing and an infinitely large number of pressure help points this expression amounts to: (57) The second, decisive force contribution results from viscosity effects (see Figure 14).Because of a lack of wall boundary conditions at the border between plugged and unplugged cells, no "zero velocity" condition can be introduced at the particle surface.What happens is that an effective "zero velocity" condition is imposed along a virtual surface including all cell centres just within the particle borders.Therefore local shear forces F Ti at the help point positions can be approximated by using the velocity value of the nearest, unplugged fluid cell u f,i , at distance h u,i perpendicular to the particle surface.Thus the overall shear force F T on the particle can be calculated as: (58) For infinitesimally fine grid spacing and an infinitely large number of pressure help points this expression amounts to: Where is the Jacobean of u f .Under the condition of incompressibility, div u f =0, this expression can be expanded and generalized.With τ being the viscous shear stress tensor, this amounts to: Pressure and shear stress contributions to the overall drag force on an ellipsoid particle are represented in Figures 13 and 14.For a limited number of discretizing surface elements N this expression yields: (62) It has to be stated that this drag-and lift force implementation is grid dependent.The applied meshes however, are structured grids with never changing resolution.An exact knowledge about particle shape and surface structure of actual dirt particles is not given.Yet plausibility commands the following statements to hold:

•
Arbitrarily shaped dirt particles rather behave non-spherically than spherically.• Dirt particles have rough surface structures than smooth surfaces.The fibre vicinity drag module gives a binary, coarse, grid spacing ∆s dependent representation of the particle surface.Consequentially the module yields significantly higher fluid-particle forces than a corresponding representation of smoothly surfaced objects.However, the qualitative drag behaviour against Re p is fine.Figure 19 summarizes the drag force behaviour of simple Stokes-flow-spheres and fibre-vicinity-module-spheres.The ratio S = ∆s/D sph is used as a parameter.
For supposedly, arbitrarily surfaced particles the CFD results can be expected to be more appropriate, than any smooth surface representation.Still, correction functions have been introduced to compensate for surface roughness, and numerical resolution effects on a user defined basis.Because of the good qualitative behaviour of the solution, the finding of a suitable correction function is comparatively simple.Possible parameters of dependence are the particle Reynolds number Re p and the grid spacing ratio S. Exemplary cases within the parameter ranges 0.05 ≤ Re p ≤ 2.0 and 0.05 ≤ S ≤ 0.5 have been evaluated.Note that the particle model is, as of now, declared valid only for creeping flow conditions Re p < 0. A consideration of Figure 17, Equ.65 and Equ.66 yields the surprising result that ζ will have to be smaller for larger S-values than for fine grid spacing, even though the shape representation gets worse.The explanation for this is given by the fact that, with larger S, the closed fluid cell volume V block decreases as compared to the analytic volume of the object V a , until S~0.5.For S-values larger than 0.5, the "large" particle model forfeits its validity anyway.
Figure 19 shows analytical, un-corrected and smoothness-corrected model results in terms of c d -values.

CONCLUSION
A non-spherical dirt particle model for filtration applications has been programmed under the Open Source CFD tool box OpenFOAM®.This work has presented the most essential, underlying concepts behind the new solver.Moreover the two drag-/lift force modules for the free flow regime and the filter fibre vicinity have been thoroughly explained and investigated.Thereby the efficient, porosity based, two-way coupling strategy was discussed.A related, sequel article [3] is to be published as well.It focuses on event force modelling in terms of particle-particle, particle-wall and particle-fibre interaction effects.In addition to that, [3] deals with data conditioning and the development of a semi-analytical result verification scheme.
Based on the geometrically versatile shape of ellipsoids, dirt particles, that feature six degrees of motional freedom, can now be modelled.A wide variety of particle shapes,anything from plates to sticks to simple spheres -, can be approximated by this concept.
The development effort behind this paper has entailed several prominent sub steps, such as: • Design of an explicit, force-and torque vector model, which reduces the modelling to the mere formulation of single force effects.

•
Introduction of an adaptive time stepping scheme for explicit Euler discretization of the PME, [4].• Device of a surface help point scheme to account for large particle effects in terms of fluid-particle, particle-fibre and particle-particle interaction.• Application of a drag force implementation that uses a combination of non-spherical, semi-empirical drag force formulas [23], a panel method to consider free flow swirling effects, a plugging method to include inter particle and particle-fibre hydrodynamics and a simple adaptation of basic concepts known from the immersed boundary method [16].• Implementation and description of an efficient particle-fluid, two-way coupling method [2].
This paper presents the framework for large, non-spherical, versatile particles with multiple abilities.The concept of additional event forces, which can be added to the PME, enables the modular addition of individual interaction effects like those, presented in [3]: The nature of the solver is such, that other, additional features could be implemented pretty easily in the future.Some possible examples are:

•
Electro static interaction • Diffusion contributions to particle motion • Particle-particle skin friction forces In a further step, the fully assembled particle module can be easily added into any laminar, continuum based, Eulerian fluid solver.
Since validation is crucial as well, extensive measures to verify solver functionality and result quality had to be taken.Article [3] reports of the development of a semi-analytical verification scheme to back up CFD results.
In the future this CFD tool will enable the virtual, purely digital pre-design of filter fibre materials.It will then be possible to upload artificially created fibre structures, conduct the CFD analysis and to use the results as a good estimate on the expected performance of the newly designed fibre material.

Figure 1
Figure 1 Digitalization by matlab utilities.Transferring stacks of grey scale images (left) to fully digitalized data (right).

Figure 4
Figure 4 Exemplary ellipsoid particle with co-rotational coordinate system.

Figure 5
Figure 5 Illustration of acting forces and torques on an assembly of non-spherical particles.

Figure 7
Figure 7 Non-spherical particle surrounded by help points and panels.

Figure 10
Figure 10 Sketch of local force balance and force effect on panel centre.

Figure 11
Figure 11 Sketch of form-and shear force contribution to panel drag force.

∑ 1 NFigure 12
Figure 12 Ellipsoid particle accelerating towards valve.Alignment along the current stream lines.Particle takes up its most stable position of least drag-and lift forces.This behaviour causes the non-spherical slip effect with relevance for filtration efficiency and particle penetration depth[3].

Figure 13
Figure13 Pressure force contribution to over all fluid-particle force.Exemplary, two-way coupled ellipsoid.Pressure build up in frontal particle area.Formation of pressure gradient across particle surface.

Figure 14
Figure 14 Shear stress contribution to over all fluid-particle force.Exemplary, twoway coupled ellipsoid.Plugging equals a zero flow velocity boundary condition at engulfed cell centres.Boundary layer is approximated.Shear stresses can be derived.

Figure 15 Figure 15
Figure 15 qualitatively shows how some two-way coupled, non-spherical particles can affect the surrounding fluid flow.To quantify the results, a fluid-particle force evaluation within the fibre vicinity module

Figure 16
Figure 16 Plot of ζ against Re p .For Re p < 0.5 there is no relevant result dependence on Re p .

Figure 17
Figure 17 Plot of ζ against S. Results (red) are fitted linearly (blue) and by a third order polynomial (green).

Figure 19
Figure 19 Plot of log(c d ) against log(Re p ). Averaged (over S range), original model results (purple) are fitted by polynomial smoothness correction, using Equ.66 (blue) to analytical Stokes drag results (yellow).Model is valid within the Stokes drag regime, log(Re p ) ≤ -0.30.