Simulation of Fluid-Structure Interaction using OpenFOAM: Filtration Processes in Deformable Media

Volltext

(1)University of Leoben, Austria Petroleum Engineering Dissertation. Simulation of FluidStructure Interaction using OpenFOAM®: Filtration Processes in Deformable Media Marianne Mataln, M.Sc. ICE Strömungsforschung GmbH, Austria &. Department of Mineral Resources and Petroleum Engineering 1st Advisor: A.o. Univ. Prof. DI Dr. techn. Wilhelm Brandstätter Department of Mineral Resources and Petroleum Engineering University of Leoben, Austria 2nd Advisor: Univ.Prof. DI Dr.mont. Thomas Antretter Institute for Mechanical Engineering University of Leoben, Austria October 2010.

(2) I declare in lieu of oath, that I wrote this thesis and performed the associated research myself, using only literature cited in this volume..

(3) Abstract [1]/id} In filtration processes it is necessary to consider the interaction of the fluid with the solid parts and the effect of particles carried in the fluid and accumulated on the solid. Traditionally, for investigation of the driving parameters, such as particle deposition and material influence, destructive tests are used. In order to provide accurate repeatable ambient requirements, simulations offer an attractive alternative. While other related publications deal with the large particle deposition model [2-5], this thesis focussed on the development of a solver to model fibre deformation effects. Pressure and traction forces, induced by fluid motion consequently lead to deformations of the solid part. According to this multi-physical problem, it is necessary to couple the differential equations of fluid motion, namely the Navier-Stokes equations and structural mechanical equations, Hooke’s law, for the solid region. For their numerical discretisation only one single computational mesh is used. This grid is changing with time and hence is recalculated at each time step to adjust to the deformation in order to conserve geometric consistency. The derived algorithm was summarized by one single solver and realised with the help of the Open Source, C++ based, computational fluid dynamics tool box OpenFOAM®. It was thoroughly validated by plausibility checks and available experimental data. With this, a strong tool for studying fluid-structure interaction phenomena on microscopic scale was developed. It was successfully applied on realistic, from CT-scans reconstructed fiber materials. Further on it was combined with the Lagrangian particle model. This provides the possibility of simultaneous simulation of all relevant physical phenomena by one single finite volume solver based on OpenFOAM®. Experiments showed a non linear behaviour of pressure drop in dependency of flow rates for soft filter materials. With the help of the newly developed filtration solver it was possible to prove this observation. Further on the particle deposition behaviour for different filter materials was investigated. New insights were gained, which underlined the high influence of the deformation of filter material on the overall filtration process. The final aim of the project is to design a filtration tool for the development and optimization of new high performance filter materials without need for performing time consuming and expensive experimental work..

(4) Kurzfassung In technischen Filtrationsprozessen sind vor allem zwei Effekte von grundlegender Bedeutung: die Wechselwirkung zwischen dem Fluid mit der Faserstruktur des Filters und die Ablagerung der, mit der Flüssigkeit transportierten Schmutzpartikel. Traditionell werden zur Untersuchung von Filtrationskenngrößen, wie der Filtereffizienz, aufwendige, destruktive Tests angewandt. Um die Zerstörung der Faserstruktur jedoch zu vermeiden und exakt reproduzierbare Bedingungen zu schaffen, bieten sich Simulationen als eine sehr gute Alternative an. Während sich andere, verwandte Publikationen mit der Entwicklung eines Schmutzpartikel-Ablagerungskonzeptes befassen [2-5], bezieht sich diese Arbeit auf die Entwicklung eines Simulationsprogrammes zur Modellierung von Deformationseffekten. Die Deformation der Filterfasern wird infolge der, auf die Oberfläche der Struktur wirkenden strömungsmechanischen Druck- und Scherkräfte ausgelöst. Im Gegenzug erfolgt eine Änderung des Strömungsfeldes durch die Bewegung des Festkörpers. Zur Modellierung dieses Verhaltens ist es erforderlich die strömungsmechanischen Differentialgleichungen des Fluides, die Navier-Stokes Gleichungen und jene der Strukturmechanik des Festkörpers, also das Hook’sche Gesetz, zu koppeln. Um die geometrische Konsistenz zu wahren, wird zur räumlichen Disketisierung ein zeitabhängiges Berechnungsgitter verwendet, welches in jedem Zeitschritt an die Deformation des Festkörpers angepasst wird. Der entwickelte Algorithmus ist in einem einzigen Simulationsprogramm auf Basis der Open Source Simulationstoolbox OpenFOAM® realisiert. Anhand von Plausibilitätskontrollen und verfügbaren, experimentellen Daten wird die entwickelte Software gründlich validiert. Auf diese Weise entsteht ein verlässliches Werkzeug zur Simulation der FluidStruktur Interaktion im mikroskopischen Bereich, welches erfolgreich an realistischen, aus CT-Scans rekonstruierten Fasermaterialen angewendet wird. Darüber hinaus wird eine Koppelung mit dem separat entwickelten Langrangen Schmutzpartikelmodell verwirklicht. Diese Verknüpfung ermöglicht somit die gleichzeitige Simulation aller relevanten physikalischen Phänomene in einem einzigen Finite Volumen Simulationsprogramm auf Basis von OpenFOAM®. In Experimenten kann ein nichtlineares Verhalten des Druckabfalles über die Dicke des Faserelements in Abhängigkeit der Durchflussraten beobachtet werden. Eine weitere, durch die Deformation beeinflusste, charakteristische Filtereigenschaft ist die Filtereffizienz, also die Fähigkeit Schmutzpartikel verschiedener Klassierungen aus dem Fluid abzuscheiden. Das neuentwickelte Simulationswerkzeug zur Filtrationsanalyse ermöglicht eine Untersuchung dieser Beobachtungen. In dieser Arbeit werden viele neuartige Erkenntnisse präsentiert, die den hohen Grad des Einflusses der Deformation des Filtermaterials auf den gesamten Filtrationsprozess unterstreichen. Die Zielsetzung dieser Arbeit ist das Design eines Filtrationswerkzeuges zur Entwicklung und Optimierung von neuen Hochleistungsfiltermaterialien ohne die Notwendigkeit von zeitraubenden und kostenintensiven Experimenten..

(5) Acknowledgement The. studies. underlying. this. thesis. have. been. performed. at. the. ICE. Strömungsforschung GmbH in cooperation with MAHLE Filtersysteme Austria GmbH. I would like to thank Prof. W. Brandstätter for giving me the opportunity to work on this thesis and the freedom to develop myself during my studies. My thanks go to all my colleagues at the ICE Strömungsforschung GmbH for their help, especially to DI Bernhard Gschaider, to whom I owe a great deal of thanks for his persistent support in every programming “case of emergency”. Special thanks are due to Dr. Gernot Boiger who was a good friend and colleague during the whole project. I have a lot of respect for his technical skills and I very much appreciated our numerous conversations. I wish to express my gratitude towards Prof. Antretter and Prof. Oberaigner from the Institute for Mechanics for the confidence they placed in me and giving me the opportunity to experience the life of a teaching assistant. Further on, my thanks go to Malcolm Werchota for his help and his talent to provide me with a lot of exciting challenges. Finally I would like to thank my family for all their support during my whole life, for their persistent help and all the wonderful time we had and still have together..

(6) Table of Contents 1. Introduction ........................................................................... 3. 2. Overview................................................................................. 8 2.1 2.2 2.3 2.4. 3. Fibre Material and Filter Types ................................................................. 8 Digital Fibre Reconstruction [2] ................................................................ 9 Particle Deposition Model ....................................................................... 11 Programming Platform: OpenFOAM® ...................................................... 13. Fluid-Structure Interaction (FSI) ........................................ 16 3.1 Introduction ......................................................................................... 16 3.2 Constraints and Specifications ................................................................ 17 3.2.1 Fluid Region ................................................................................... 18 3.2.2 Solid Region ................................................................................... 18 3.3 Governing Equations ............................................................................. 19 3.3.1 Governing Equations for Fluid Flow ................................................... 19 3.3.2 Governing Equations for Solid Deformation ....................................... 21 3.4 A Single Algorithm................................................................................. 23 3.5 Discretisation ........................................................................................ 26 3.5.1 Spatial Discretisation of the Computational Domain............................ 27 3.5.2 Equation Discretisation .................................................................... 28 3.5.3 Moving the Computational Mesh ....................................................... 29 3.6 Structure of an OpenFOAM® Case for FSI Simulation ............................... 33 3.7 Application to a Simple Filter Model ........................................................ 38 3.7.1 Boundary Conditions ....................................................................... 39 3.7.2 Results ........................................................................................... 40 3.8 Further Improvement ............................................................................ 44 3.8.1 New Boundary Conditions ................................................................ 44 3.8.2 Collision Algorithm .......................................................................... 51 3.8.3 Combination ................................................................................... 56 3.8.4 Investigation of Different Branching Types ........................................ 63. 4. Application to Digitally Reconstructed Fibres ..................... 65 4.1 Computational Optimisation ................................................................... 66 4.1.1 Parallelisation ................................................................................. 66 4.1.2 Restart ........................................................................................... 71 4.2 Results ................................................................................................. 73 4.2.1 Fibre Deformation ........................................................................... 73 4.2.2 Combination with Dirt Particles......................................................... 79. 1.

(7) 5. Influence of FSI on Overall Filter Efficiency ....................... 83 5.1 Experimental Setup ............................................................................... 85 5.2 Comparison of Results of Experiments and Simulation .............................. 87 5.2.1 Nonlinear Pressure Drop Behaviour................................................... 87 5.2.2 Filter Fibre Efficiency Curves ............................................................ 90. 6. The Graphical User Interface [2] ......................................... 94. 7. Conclusions and Outlook ..................................................... 97. 8. Nomenclature ....................................................................... 99. 9. List of Figures ..................................................................... 101. 10. References .......................................................................... 107. 2.

(8) 1. Introduction. Filtration processes are important and used in a wide variety of industries. Examples therefore include the automotive industry, sewage filter systems, or even facilities present in our household, such as the coffee machine. This thesis will deal with oil filters, which are in widespread use within the automotive industry, Figure 1.1. Motor oil is used to lubricate the components of the engine. Additional functions include reducing wear, cooling the engine and protecting it from corrosion. Acting by itself, the oil would soon become saturated with dirt particles, such as dust that is ingested into the engine, as well as abrasive metal parts. This will consequently lead to wear of internal parts, reduced protection from corrosion and in the worst case to engine failure. At this point oil filters come into play, as only clean oil and proper flow can guarantee continuous excellent engine performance. Oil filters consist of a strong steel box, that can withstand high oil pressure (up to 6 bar), an anti-drain back valve, that shall create minimum backpressure, a pressure relief valve, that doesn’t leak and a strong filter media element. Microscopically seen, this media consists of a maze of different fibres, which will capture most types of dirt particles. The surface of the filter media is folded to extend its size in order to trap a substantial amount of particles and to stabilize the overall filter media. As time goes by, dirt particles will clog the filter, restricting oil flow and hence leading to oil starvation. As a consequence less supply of oil will damage and in the worst case destroy the engine. In order to extend lifetime and reach higher efficiencies of filters, it is sensible to improve these in every possible form. Generally these improvements require experimental runs, which tend to be time consuming and very costly. In the majority of cases those experiments are carried out with destructive tests. In that case filters are cut or burned to gain knowledge about the investigated quantities. Yet these tests cannot be repeated in order to obtain other characteristics of the filter or to prove a different outcome. Another disadvantage of experiments is the difficulty to investigate filter material in detail, i.e. to find out, which influence the deposition of particles or the deformation of fibres would have had on the overall filter characteristics. Therefore CFD simulations can offer an attractive alternative to experiments. Apart from the cost saving fact, every quantity of a filter can 3.

(9) independently be investigated, displayed and processed for further usage. The virtual tests can be repeated as often as necessary with exactly the same initial situation and properties. CFD simulations on microscopic scale enable the creation and testing of new filter materials without extensive experimental runs.. Figure 1.1: Detailed view of an oil filter In modelling of such a filter material, first it is required to investigate the main governing factors of the filtration process. One of those is that the flow of the fluid through this filter induces a pressure distribution on the fibres, which leads to a deformation of the solid parts. This deformation results in a different assembly of the fibres, which will have significant impact on the other factors, as well as on the overall material. Another effect results from the fact that the fluid is carrying abrasive particles resulting from engine wear. The efficiency is measured on the ability of the filter material to absorb different particle sizes. It is essential to understand the underlying driving factors for further improvement of filter materials. Due to the complexity of the issue, a change in fibre morphology (i.e. pore size diameter) cannot be linearly linked to i.e. filter efficiency, because it influences the whole hydrodynamic situation. This thesis represents an extensive attempt to create a tool which can increase the understanding of filter effects and dynamic parameter dependencies by means of computational engineering and simulation technology. A detailed, deterministic calculation model, which simulates the most important filtration effects on a microscopic level, has been created. Figure 1.2 sketches out the basic concept behind this novel scheme. [2]. 4.

(10) In a firrst step, computer-t c omographyy (CT) sca ans are co onducted o on “real liffe” filter fibre sa amples. Th he CT output data is compiled in stacks off two dime ensional (2 2D) gray scale im mages of the fibre. Then T the da ata is read in, digitaliized, and p processed to a full 3D reconstruction n of the microscopic m c filter elem ment. The e 3D objecct is autom matically d by a stru uctured grrid generattor, in ord der for the geometryy to be utilized as meshed bounda ary framew work for on ncoming CFD C calcula ations. Thiss is where e the resultt of the main developmen d nt task com mes into play: p A CF FD tool, de esigned and program mmed in order to t resolve the dynam mic filtratio on situation for a usser definab ble set of process variable es, within the recon nstructed fibre f element. Produ uced simulation resu ults can then be e used to estimate e th he perform mance and suitability s o the tested medium of m. [2]. Figure 1.2: Sketcch of the principle p simulation concept. CT C scans yyield stacks of 2D grey sccale image es (left), which w are transferred t d into 3D reconstrucctions of th he fibre (middle e). The 3D images arre meshed and provide the geo ometry for the CFD filtration f solver to t be creatted (right)... Constant checks for result plausibilityy and validation have e to be inte egral partss of any elopment effort. e In order o to qu ualitativelyy and quan ntitatively validate v seriouss CFD deve the ressults, an exxtensive exxperimenta al set up has h been crreated and d a semi-empirical validatiion scheme e has been n devised. Figure 1.3 3 gives an overview of the und derlying develop pment sch heme whicch links th he experim mental- and d the sim mulation sid de. The application of thiss method leads l to a continuous adjustme ent and im mprovemen nt of the CFD mo odel, accorrding to the equivale ent, experim mental resu ults. [2]. 5.

(11) Figure 1.3: Experrimental- and CFD de evelopmentt From the t beginn ning the de evelopmen nt project was divide ed into four major working w areas, as a seen in Figure 1.4 4 [2] ( Digital Fibre Recconstruction (DFR) from (1) f CT in nformation n, as well as the suitable e meshing of the 3D data. ( Creatio (2) on of a Flu uid Structure Interacttion (FSI) tool in ord der to han ndle the fibre de eformation n effects un nder the influence of fluid flow. ( Development of a detailed,, dirt particcle- and deposition m (3) model, cap pable of simulatting spherrical- and non-spherical dirt particle behaviour in- and outside e of the miccro scale fibre vicinity y. ( Validatiion of simu (4) ulation ressults. Devisse of an ap ppropriate,, experime ental set up to verify v solve er function nality and to t provide additionall insight in nto filter fibre be ehaviour and charactteristics.. 6.

(12) Figure 1.4: Overview of th he four ma ajor areas of develo opment behind the filtration f T development of a suitable e fluid-stru ucture inte eraction so olver for solver project. The fibre de eformation n is at the focus f of this thesis. , 2004 29 /id} o rela ated publiccations [2-5] extensively deal with the d developme ent of a While other, detailed d, dirt partticle- and deposition n model an nd the DFR R utility, th his thesis will w only briefly discuss those two su ubjects. Th his thesis will be foccussing on the development of the fluid-struccture intera action parrt. The fully optimize ed and im mplemented d solver constitu utes a stro ong tool fo or the inve estigations on how the t charactteristics off the oil filter will w change e with the deformation of the fibres due e to oil flo ow. Furthe er on in chapterr 5 the con nducted ve erification process p willl be discusssed and d displayed in n detail.. 7.

(13) 2. Overv view. 2.1. Fibre Material M l and Fillter Typ pes. Early designs d of oil o filters co omprise stteel wool, wire w meshes and me etal screen ns. Over the yea ars disposal filters have becom me popula ar and with h that, cellulose and papers were used. u Cellu ulose is a natural ma aterial and d presents a random m irregular field of fibres to o the oil.. Figure 2.1: Detailled view off fibres from a filter media m In gene eral, engin ne oil filterr material has h to yield the follo owing feattures to maximize service e intervals: •. ossible surfface to trap p as many particles as a possible e Largest po. •. Least posssible restricction to oil flow. •. Handle as much conttamination n loading as possible,, before ne eed to be replaced r. •. er not to afffect perforrmance of the vehicle e Light weight, in orde. t ideal balance of dirt d trappin ng efficienccy and dirtt holding ca apacity. It is vittal to find the Only th hrough thiss, excellentt engine prrotection will w be guarranteed. For be etter qualityy filters, cellulose c iss blended with cotto on or micrroscopic sy ynthetic fibres. Synthetic media havve smaller passages that trap more partiicles, none etheless o the oil flow, since the fibrres are thinner. Figure 2.2 showing less restriction to illustrattes the diffference.. 8.

(14) Figure 2.2: Difference between cellulose and synthetic filter media Another way of enhancing filter material performance is to use the thinnest possible fibres. As a matter of fact, the thinner the fibres, made from a predefined amount of material are, the longer these will be and hence possessing a larger effective surface area. They can trap a larger amount of particles and additionally show much less restriction to the oil flow. This can also be realised by using synthetic fibres. The advantage of those is the possibility of specific manufacturing of the fibres regarding shape and thickness. A special type of filter is the so-called “depth” filter. It has a passage size gradient, which means that it traps different size of particles at different spots in the media. The deeper the element, the smaller the passages will be and the smaller the particles which are trapped. Its main advantage is that it will hold more particles before blocking and needing to be replaced. High end oil filters use fibre glass or extremely fine metal mesh. This will enhance the stability of the media.. 2.2. Digital Fibre Reconstruction [2]. The ability to realistically model micro scale filtration processes in filter fibre materials is in large part based upon the realistic reconstruction of micro scale filter fibre geometries. Within the context of the development effort behind this work, a sophisticated method to digitally recreate real geometries was applied. In a first step, computer tomography (CT) scans are conducted on the fibre material to be. 9.

(15) investigated. The data yielded by the CT scans are stacks of 2D grey scale images seen in Figure 2.3 (left). MatLab® [6] based reconstruction algorithms have been programmed in order to process the CT data. The picture stacks can be uploaded and the individual slices are then analyzed. Local picture areas of higher grey scale intensities are recognized as fibre regions, which can be clearly distinguished against the low-intensity background. Identified fibre slices are then quantified, their pixel area is calculated and their local centres and radii are determined. By applying a skeleton algorithm [7] the centres of consecutive pictures are interconnected to constitute the basic, local fibre framework. By applying the calculated radius information attached to each centre point, the actual fibre structure is recreated as a 3D digital data matrix. This can be visualized as seen in Figure 2.3 (right).. Figure 2.3: Fibre reconstruction and digitalization by MatLab® utilities. Stacks of grey scale images (left) out of CT scans are transferred to fully digitalized data matrices (right). In a next step the digital data is automatically discretised into a structured, hexahedral grid mesh with a user definable cell-spacing-to-pixel ratio. This means that, if the CT scan resolution can be kept constant, a uniform spatial discretisation rate for any filter fibre simulation can be guaranteed. Thereby one of two modes of spatial resolution can be chosen: Either the finer mode, which features a spatial resolution of 1.6µm/Pixel or the coarser mode, which features a resolution of 3.2µm/Pixel. The reconstruction utility yields perfectly interfaced grids, of both the 10.

(16) fluid an nd the sollid region of the fib bre sample e. Figure 2.4 2 presen nts an exe emplary, structured, micro scale fibre e grid mesh.. Figure 2.4: Filter fibre samp ple discretiised into a structured d fluid- and d solid hexxahedral grid me esh. Dimen nsions: 200 0µm x 200µm x 300µ µm. Numbe er of cells: ~6.0 x 10 05 epared, strructured grid meshess serve as geometry boundary conditionss for the The pre simulattor to be developed.. 2.3. Particlle Depossition Model M. The se econd part of the prroject wass to develo op a tool for simula ation of large dirt particle es. Low de ensity clou uds of dirrt particless, ranging from 5 tto 50 miccrons in diamete er, occur in the oil stream. They T are classified as “large particles””, which means that their size is much larger than t the size of a grid cell in th he simulatio on grid. Those large particles conta ain a certa ain mass and a momentum and interact with w the flow fie eld. This is in contrrast to sim mulation off “small pa articles”, w where usually the temporral and spa atial evoluttion via the solution of Euleria an conserva ation equa ations is sought. The difference is sh hown in Fig gure 2.5.. 11.

(17) Figure 2.5: Particcle simulation with sm mall particle es (left) an nd large pa articles (rig ght) The ma ain effectss governing particle motion arre wall intteractions, interactions with filter material m an nd the colllision betw ween them mselves and the plug gging effe ect. The latter effect e imp plies that with time e, more and a more particles find them mselves entangled in the e fibre stru ucture and d slowly acccumulate there, thus increassing the effectivve solid su urface and causing plugging p off the flow. Even though those effects ic filter cha occur at a a microsscopic leve el, they will change macroscop m aracteristiccs, such as poro osity, presssure drop p and ove erall filter efficiency.. In orderr to simulate the encoun ntered phenomena, not n only is the implem mentation of a detailed particle e model necessa ary, but further f the e calculatiion of particle force interacttions is esssential. Therefo ore a discrrete phase Lagrange model wass develope ed, which ttakes into account a the two o-way coup pling betwe een the flu uid and the e particles. The particcles are related to Lagrang ge objectss. This imp plies that they are regarded as individu ual objectss. Their behavio our, such as velociity, accele eration and position n, is totally depend dent on surroun nding fluid flow conditions such h as flow velocity and d flow presssure.. Figure 2.6: Particcle cloud in n a digitallly reconstrructed fibre e geometryy. Large numbers n e impact events e occu ur. [2] of partiicle-particle 12.

(18) The mo odel also enables e the simulatio on of vario ous shapess of dirt pa articles. Th hey can be of regular sp pherical shape, sho own in Fig gure 2.6, and also of nonsp pherical, ure 2.7. For F developing a rea alistic dirt particle m model it is vital to presentted in Figu regard their morp phologies. For examp ple, a sphe erical particcle model w would significantly undere estimate flu uid skin fricction and form f drag forces. f Oth her relevan nt implicatiions are that elo ongated, non-spheric n cal particle es show the tendencyy to align themselves in the fluid stream thus increasing g the possibility of slipping thro ough a porre, wherea as discs, gh rotation nal relaxatiion times, more likely y block porres. with hig. Figure 2.7: Nonsp pherical dirt particless entering a realistic digitally re econstructe ed fibre structure [2] If those e and man ny addition nal effects are not ta aken into account, a th he pressurre drop, develop ping after some time e, may be e over- or underestim mated and d therefore e only a nonsph herical partticle modell will allow a realistic filter chara acterisation n. d de escriptionss regarding g the bacckground and devellopment of o large Very detailed sphericcal and non nspherical particle models for filtration f sim mulations are laid ou ut in [25]/id}. 2.4. Progra amming g Platforrm: Open nFOAM®. The en ntire CFD related so oftware de evelopmen nt behind this thesis was con nducted within the framew work of th he Open Source CFD D package OpenFOAM M® (Open Source Field Operation O a Manipu and ulation) [8--10]. It is a CFD toolbox, which is applicable for a wide e range of o problem ms in conttinuum me echanics. It provide es the user with 13.

(19) standard solvers, utilities and libraries, whereas libraries are repositories of function related software tools that can be accessed by solvers and utilities. All can be specifically selected according to the governing physics of the problem. The underlying programming language is the object oriented language C++ [11,12]. The syntax is similar to the notation of the differential equations being solved. For example, the Navier-Stokes equations for incompressible flow read:. ∂(ρui ) ∂(ρui u j ) ∂p ∂ 2 (η ⋅ ui ) + =− + 2 ∂t ∂x j ∂xi ∂x j. 2.1. Where ui and uj are the velocities along the coordinate directions xi and xj, ρ is the density, p is the pressure and η is the dynamic viscosity. [13] This is implemented into the source code as following: solve (. fvm::ddt(rho, U) + fvm::div(phi, U) - fvm::laplacian(eta, U) == - fvc::grad(p). );. where φ = ρ ⋅ U (written as “phi” in the code) The source code of the program has been made Open Source and thus is publicly available to anyone under the constraints of the GPL [14]. Therefore every aspect of the underlying source code can be altered as required by the user, which even permits the creation of whole new solvers, if necessary. It allows easy and direct implementation of new software modules at any point in the program. It is easy to automatise with the help of scripts and it extensively uses generic features (Templates). These several advantages make OpenFOAM® highly efficient and a very flexible tool, which in addition is free of license costs. The only disadvantage of OpenFOAM® is that it is not provided with a Graphical User Interface. This means that all inputs have to be provided by means of text files, which requires a higher effort for familiarising oneself with the software and a significant amount of prior knowledge of physics and programming techniques. Nevertheless OpenFOAM® is an 14.

(20) appropriate environment to create new solvers and hence model multiphysics problems such as the fluid-structure interaction found in oil filter applications.. Pre-and Postprocessing: In this thesis the meshing was conducted via the commercial FLUENT® mesh generator GAMBIT® [15] or via self written meshing utilities. All results were post processed and visualized with the Open Source visualization tool ParaView [16] by Kitware®.. 15.

(21) 3. FluidStructure Interaction (FSI). 3.1. Introduction. In modelling of filter media a very important physical effect has to be considered concerning the interaction between the fluid and the fibre material. The fluid flow induces forces on the solid regions, which lead to deformation of the fibres and hence to changes of the overall material structure under working conditions. This induces severe modifications of the permeability of filters. The main challenge now is to simulate this situation with all the underlying physical phenomena considered. Traditionally, for modelling fluid flow, CFD (computational fluid dynamics) codes are widely used, whereas for mechanical stresses CASA (computer aided stress analysis) codes are commonly employed. In the majority of cases, CFD codes are based on the so-called Finite Volume (FV) Method, where on the other hand CASA codes are mostly based on Finite Element (FE) principles. This procedure requires that two independent computer codes be used to calculate both, fluid flow and solid-stress simulation for the same system. Furthermore a third code is required for coupling and data management. This implies that all codes must be available and the user must have knowledge of all three of them. Another problem is that the CFD code must be converted into an appropriate format readable by a CASA program and vice versa. For both circumstances this conversion will cause a certain degree of loss in accuracy, especially as the two-way information exchange has to be undertaken iteratively, until convergence is achieved. The solution is to combine the overall calculations in one single computer code. To enable the coupling, similarities between the fluid flow equations and solid stress equations have to be found. This will give the possibility to devise an algorithm, which will solve the solid-stress equations in one part of the field and the fluid-flow ones in another, yet realise both in one single computer code. The question now is which basic principle shall be used, the Finite Element (FE) method or the Finite Volume (FV) method.. 16.

(22) The FE method provides an accurate and stable result for transient simulations of structural deformation. It is based on the variational principle [17] and uses predefined shape functions dependent on the topology of the element. It easily extends to higher order discretisation, produces large block matrices with usually high condition numbers and as a consequence relies on a direct solver. The discetisation is non conservative and hence does not guarantee to satisfy the conservation equations of the fluid flow. In case of non linear equations or discontinuous coefficients the FE method may show instabilities. Details can be found in [18,19] The FV method can be considered to be a particular case of the FE method. For the FV method, the shape functions are regarded as piecewise linear and hence allows a conservative discretisation. Therefore it can handle complicated, coupled and nonlinear differential equations, widely used in fluid flows [20]. The non-linearity is treated in an iterative way and creates diagonally dominant matrices well suited for iterative solvers. From this follows that the FV method also allows discontinuous solutions, such as they appear in local mesh refinement with discontinuous mesh line intersections. Additionally it has minor demands on the quality of the computational mesh. The decisive factor was the available environment for the development of a new solver. Most of the conventional simulation programs do not allow access to the underlying computer code. OpenFOAM® being an Open Source program permits a direct access to every part and equation implemented in the program and additionally allows for modification. Due to these advantages, OpenFOAM® is the best suited environment for development of a new solver for fluid-structure interaction. It is based on the FV method, which additionally shows the best behaviour for simulation of filtration processes.. 3.2. Constraints and Specifications. For simulation of oil filters not all factors, which govern fluid flow and solid deformation, have to be considered. A few simplifications can be made in order to reduce the time and effort required for CFD simulations. 17.

(23) 3.2.1. Fluid Re egion. The sim mulation iss carried ou ut at a miccroscopic scale s with the intenttion of stattistically averaging resultss for several regions in order to t retrieve macroscopic information in later ph hases of developmen d nt. Due to the size of o the geom metry, which is in th he scale of a few w microns,, the Reyn nolds numb bers are ve ery low and d range be etween 0.5 5 and 5. In that case the flow f can be classified d as creepiing flow an nd turbulen nce effectss do not t be takken into account. a T Temperatu re also has h little influence on the have to deform mation of th he filter material and d the existting particles and the erefore the e model can be assumed to be isoth hermal. Th here is only y oil to be considered d as a sing gle fluid ence single e phase flo ow assump ptions are valid as well. w Genera ally not ve ery high and he pressurres are occcurring du uring the filtration process p an nd thus ch hanges of the oil densityy are not significant. s . From thiss follows, that the flow f can b be assumed d to be incomp pressible.. 3.2.2. Solid Re egion. The efffect of com mpression of the ovverall filter material is i more im mportant th han the deflectiion of a sin ngle fibre. Hence it is satisfactory to rega ard only sm mall deform mations of fibre es, where the longittudinal mo ovement u is much less than the deflecction w. Therefo ore u can be b neglecte ed.. Figure 3.1: Deforrmation of a single fib bre due to an arbitrary force F It that case, the deformation behaviiour is line ear, which means th hat the cha ange of load off a certain factor ressults in a deformatio on of the same facttor. Thereffore the stresse es are pro oportional to strainss accordin ng to Hoo oke’s law. Additiona ally the materia al can be considered c d as elasticc and hencce, there only o are d deflections and no plastic deformatio ons of the e fibre itse elf. Finally no initial stresses w will be tak ken into account. 18.

(24) 3.3. Governing Equations. 3.3.1. Governing Equations for Fluid Flow. The basis for deriving equations for fluid flow is mass and momentum conservation. Newton’s second law states that the rate of change of momentum of a fluid particle is equal to the sum of all forces acting on the particle. All of the following equations are presented in Einstein’s notation. The rate of change of momentum in a unit cell is given by:. ρ. dui ∂u ∂u = ρ i + ρu j i ∂t dt ∂x j. where ui and uj are the fluid flow in coordinate directions xi and ρ the density. called substantial derivative,. 3.1 du i is dt. ∂u ∂u i is the partial derivative and u j i is the convective ∂t ∂x j. term [13,21-23]. There are two types of forces acting on a differential fluid unit cell: •. Surface forces: [21] ¾ Pressure forces (the negative sign means it is pointing outwards of the cell volume):. −. ∂p ∂xi. 3.2. ¾ Viscous forces evolve due to stresses applied to the control fluid volume:. ∂τ ij ∂x j •. 3.3. Body forces: ¾ Source term:. ρg i where p represents the pressure on the fluid cell, gi the gravity vector and. 3.4. τ ij the. normal and shear stresses. 19.

(25) Adding up these contributions, we derive the following: 2005 10 /id}. ρ. ∂ui ∂u ∂p ∂τ ij + ρu j i = − + + ρg i ∂t ∂x j ∂xi ∂x j. 3.5. For Newtonian fluids it is valid that there is a linear relationship between stress and strain rate. The proportionality factor is constant at moderate temperatures and is called the dynamic viscosity η of the fluid. For one dimensional flow one obtains:. τ =η ⋅. ∂ui ∂x j. 3.6. Stokes extended Newton’s idea from a simple 1D flow to a multidimensional flow leading to the Stokes’ relations:. ⎛ ∂ui. τ ij = η ⋅ ⎜⎜. ⎝ ∂x j. +. ∂u j ∂xi. −. 2 ∂u k ⎞⎟ δ ij ⎟ 3 ∂xk ⎠. 3.7. where δij is called the Kronecker symbol and is defined as. ⎧1 ⎩0. δ ij = ⎨. für i = j für i ≠ j. 3.8. By inserting Eqn. 3.7 into 3.5 we derive. ρ. du i ∂p ∂ =− + ∂xi ∂x j dt. ⎛ ⎛ ∂u ∂u j 2 ∂u ⎞ ⎞ k ⎟⎟ ⋅ ⎜η ⋅ ⎜ i + − δ + ρg i ⎜ ⎜ ∂x j ∂xi 3 ij ∂xk ⎟ ⎟ ⎠⎠ ⎝ ⎝. 3.9. After rearranging the second term on the right hand side by using Schwarz’ theorem regarding the permutability of derivatives [24] and applying Eqn. 3.1 we derive the fundamental Navier-Stokes equations for conservation of momentum:. ∂ui ∂ui ∂ 2ui ∂p ⎛ 2 ⎞ ∂u k ρ η + ρg i + ρu j =− +η ⋅ + ⋅ ⎜1 − δ ij ⎟ ⋅ 2 ∂t ∂x j ∂xi ∂x j ⎝ 3 ⎠ ∂xk. 3.10. In oil filtration processes the flow is regarded as incompressible, which means that δuk/δxk = 0. Further on, due to simulation on microscopic scale, the influence of the gravity is negligible.. 20.

(26) Hence, the Navier-Stokes equations for incompressible flow read as follows:. ∂u i ∂u i ∂ 2ui ∂p + ρu j =− +η ⋅ ρ 2 ∂t ∂x j ∂xi ∂x j 3.3.2. 3.11. Governing Equations for Solid Deformation. The mathematical model describing small deformation of solids is based on the three dimensional stress distributions. From Newton’s balance equation of momentum follows that. ∂ 2 ui ∂ (σ ij ) + X bi ρ 2 = ∂x j ∂t. 3.12. where for i = j σij represent the normal stresses, for i ≠ j σij = τij represent the shear stresses, ui the displacement vector and Xbi are body forces. This is based on the Cauchy theory, which says that the state of stress at a certain point in a body is completely defined by the nine components (six independent) σij of a symmetric second order Cartesian tensor called the Cauchy stress tensor.. ⎛ σ xx τ xy τ xz ⎞ ⎟ ⎜ σ = ⎜ τ yx σ yy τ yz ⎟ ⎟ ⎜τ ⎝ zx τ zy σ zz ⎠ Figure 3.2: Three dimensional stresses on a solid element (left) with the Cauchy stress tensor (right) For elastic deformation the relation between stress and strain in case of an isotropic homogenous structure is linear.. 21.

(27) Figure 3.3: Relation betwee en stress and strain acting a on a solid bodyy For the e elastic do omain the generalized g d Hooke’s law is valid d:. σ ij = E ijkl ⋅ ε kl. 3.13. where σij is the stress disttribution, Eijkl is the elasticity tensor t and d εkl is the e elastic strain distribution d n. In the special s case of an iso otropic ma aterial, the symmetry y of the fourth order elastticity tenso or can be taken t adva antage of. The tenso or can be reduced r to be dependent d on only tw wo constants. In component no otation the constitutivve law can be written n as. λ λ ⎡σ 11 ⎤ ⎡2μ + λ ⎢σ ⎥ ⎢ 2μ + λ λ 2 ⎥ ⎢ 22 ⎢ λ ⎥ ⎢ λ ⎢σ 33 λ 2μ + λ 3 = ⎥ ⎢ ⎢ ⎢σ 12 ⎥ ⎢ ⎢σ ⎥ ⎢ 0 ⎢ 13 ⎥ ⎢ ⎢⎣σ 23 2 ⎥ ⎦ ⎢⎣. 0 2μ 2μ. ⎤ ⎡ε 11 ⎤ ⎥ ⎢ε ⎥ ⎥ ⎢ 22 ⎥ ⎥ ⎢ε 33 ⎥ ⎥⎢ ⎥ ⎥ ⎢ε 12 ⎥ ⎥ ⎢ε 13 ⎥ ⎥⎢ ⎥ 2μ ⎥⎦ ⎢⎣ε 23 ⎥⎦. 3.14. wherea as εij is defined by. ε ij =. 1 ⎛⎜ ∂ui ∂u j + ⋅ 2 ⎜⎝ ∂x j ∂xi. ⎞ ⎟ ⎟ ⎠. 3.15. 22.

(28) λ and μ are called the Lamé constants and have the following relationship. E 2 ⋅ (1 + ν ). 3.16. ν ⋅E (1 + ν )(1 − 2ν ). 3.17. μ=. λ=. where E is the Young’s modulus and ν is called the Poisson number. Additionally μ = G is called the shear modulus. After combining the equations (3.13 - 3.17), we derive Lamé-Navier’s equation for solid displacement: 2 ∂ 2ui ∂ 2ui 2ν ⎞ ∂ u j ⎛ ρ 2 =G⋅ + G ⋅ ⎜1 + + X bi ⎟ 2 ∂t ∂x j ⎝ 1 − 2ν ⎠ ∂xi ∂x j. 3.4. 3.18. A Single Algorithm. The interfacing between two different codes is rendered unnecessary, when it is possible to develop a single computer code, which can then solve the solid-stress equations and displacements in one part of the field and other one for fluid-flow, i.e. fluid velocity in another. If a closer look is taken at the governing equations for solidstress and fluid flow, derived in chapter 3.3, similarities are detected. Therefore it is possible to couple the two equation systems (Eqn. 3.11 and Eqn. 3.18) and develop one single computer code, as shown below. The main difference is the pressure gradient term of the Navier-Stokes equations, Eqn. 3.11, which is absent in the displacement equations for solids. This term can be seen as a body force from external sources represented by the term Xb in the solid displacement equations.. 23.

(29) Genera al Navier Sttokes’ Equa ations (3.1 10):. ∂ 2 ui dui ∂p ⎛ 2 ⎞ ∂u k +η ⋅ ⋅ + ρ η =− ⎜1 − δ ij ⎟ ⋅ 2 ∂xi dt ∂x j ⎝ 3 ⎠ ∂xk Lamé-N Navier’s Eq quations (3 3.18): 2 ∂ 2ui ∂ 2ui 2ν ⎞ ∂ u j ⎛ ρ 2 =G⋅ + G ⋅ ⎜1 + + X bi ⎟ 2 ∂t ∂x j ⎝ 1 − 2ν ⎠ ∂xi ∂x j. This wa ay a coupling betwee en the fluid d flow and d solid-stre ess equatio on systemss can be directlyy achieved.. The presssure field of o the fluid region can be mapp ped as a bo oundary conditio on on the solid regio on and hen nce added as an exte ernal force to the terrm Xb in Eqn. 3.17. It indu uces deform mations in the solid regions, r i.e e. the fibre material. The forces acting on the fibre mate erial can be decomposed into o pressure e based (norma al) and a sh hear stresss based (ta angential) forces, f as shown in F Figure 3.4.. Figure 3.4: Sketch of the fo orces acting g on a surfface by the e flow of a fluid In tenssor notatio on, the pre essure forcce Fpi actin ng on an interface A given its normal vector ni is define ed as. F pi = ∫ p ⋅ ni ⋅ dA. 3.19. A. In addiition to the e pressure force Fpi th he flow of the fluid with w the fre ee stream velocity u∞ also o induces shear forcces Fτi on the solid region. These T force es are cau used by velocityy gradientss normal to o the solid surface an nd also have a major influence e on the deform mation of a solid:. 24.

(30) Fτi = −η ⋅ ∫ ∂ i u j n j dA. 3.20. A. where η is the dyynamic visccosity and ∂ i u j is the e velocity gradient g fie eld. It is apparent that ussing only th he normal pressure field f as a boundary b co ondition on n the deformation of the solid s part is not sufficcient to achieve physsically plau usible results. Hence the two ve ectors can be summa arized and written in differentia al form, wh hich is. ∂ (Fpi + Fτi ) = p ⋅ ni − η ⋅ ∂ iu j n j ∂A. 3.21. ~ ⎛ d (Fpi + Fτi )⎞⎟ = (∂ i p − ηΔui ) ⋅ ni X b = ∂i ⎜ ⎠ ⎝ dA. 3.22. Furtherrmore:. From th his we derive the sou urce term, which is in n tensor no otation. X bkk =. ∂ ~ X b = ∂ k p − η ⋅ Δu k ∂nk. 3.23. All equ uations and informa ation can be summa arized in one single e algorithm m. It is displayed in Figurre 3.5.. Figure 3.5: Algoriithm to mo odel Fluid/S Structure interaction. 25.

(31) There are a two re egions defined, one for f the fluiid and one e for the sttructure. First F the Navier--Stokes eq quations arre solved in i the fluid d region. From F this the pressu ure and velocityy field is derived. d At the interfface betwe een the sollid and the e fluid region, the pressurre and tracction incre ement (δp and δt) arre transferrred as boundary conditions from th he fluid side to the so olid side (seen in Figu ure 3.6) an nd the corrresponding g source term is calculated d as derive ed in Eqn. 3.23. 3. Figure 3.6: Transsfer of coup pling data In the next step p, the Lam mé-Navier equations e are solved d for the solid regio on. The resultin ng displace ement incre ement (δuD) of the solid is passsed back to o the fluid region. The de erived algorrithm is im mplemented d and integ grated in an a explicit way. Essential for the succcessful ap pplication of o the algo orithm is an efficientt method tto compenssate for the de eforming geometry g o the sollid region,, i.e. a moving of m com mputationa al mesh strategy for the fluid f region n. This will pass the informatio on of the ssolid displa acement back to o the fluid region. This topic willl be discusssed in the e following chapter.. 3.5. Discrettisation n. There are a hardlyy any analyytical solutions for the t derived d equation ns in chapter 3.4. Therefo ore those equations have to be b solved with w numerical methods. In this case, the Fin nite Volume e method is applied.. The overrall domain n of intere est is decom mposed into a finite f numb ber of conttrol volume es (CV). In n each cell the conserrvation equ uations, as theyy are for mass, m momentum and d energy are valid. The T purposse of discre etisation is to tra ansform th he partial differential d equationss into a corrrespondin ng set of algebraic equatio ons, which can be solved easilyy [25-27].. 26.

(32) The procedure is separated into two parts: spatial discretisation of the computational domain and equation discretisation.21]. 3.5.1. Spatial Discretisation of the Computational Domain. This part comprises discretisation of time and space. For space discretisation the overall spatial domain is subdivided into a number of discrete polyhedral control volumes (CV), which fill the overall volume and do not overlap. Each CV is bounded by a set of flat faces and each face has only one neighbouring CV. The cell faces of the mesh are divided into internal faces, which delimit one CV from the neighbouring one and boundary faces, which constitute the boundary field of the overall domain. In contrast to the Finite Element Method, in finite volume techniques the topology of the volume elements is not important as there are no topology dependent shape functions necessary. Hence different cells of general polyhedral shapes with a variable number of neighbours can be used, as displayed in Figure 3.7.. Figure 3.7: Typical 3D finite volume elements This creates an arbitrarily unstructured mesh on which the governing equations are subsequently solved. For time discretisation the time interval is split into a finite number of time-steps. This means that the time derivative d/dt of the conservation equations is discretised with the help of those discrete time steps. For more detailed information on those topics see [26-28]. 27.

(33) 3.5.2. Equation Discretisation. The equation discretisation produces a numerical description of the computational domain. In FV methods the basis for discretisation is the integral form of the equations, i.e. the Navier-Stokes’ equations for momentum conservation (Eqn. 3.11):. ∂ui ∂ (ρui ) ∂ ρ dV + ⋅ u ∫CV ∂t ∫CV j ∂x j dV = CV∫ ∂x j. ⎛ ∂ui ⎜η ⋅ ⎜ ∂x j ⎝. ⎞ ⎟dV + Sφ dV ∫ ⎟ CV ⎠. 3.24. where S φ represents an arbitrary source term. Gauss divergence theorem claims that changes inside a certain volume are equal to the fluxes over the boundary surfaces:. ∂ui ∫CV ∂x j dV = ∫A ui ⋅ ni ⋅ dA. 3.25. From this, the change of a quantity can be calculated by considering the fluxes over the boundary surface of the CV. Applying Eqn. 3.25 we obtain the integral form of the transport equation:. ∂ui ∂ ( ) ρ + ρ ⋅ ⋅ = η ⋅ u dV u u n dA ∫ i ∫A i j i ∫A ∂x j ni dA + CV∫ Sφ dV ∂t CV. 3.26. From this point on the equations are solved in a segregated manner. That means that each component is treated separately. After discretisation of the equations a system of linear algebraic equations is obtained. It describes the change of a system over time. This set of equations has to be solved in an iterative manner. After each iteration step, the linear equation system is calculated again on the basis of the previous iteration. This will be done until a convergence criterion is reached. A very important constraint is that the algebraic equations are solved in a way that the overall integral balances, i.e. the mass and momentum conservation, are fulfilled. This is called “conservative” discretisation. For more detailed information on discretisation of equations see [26,27]. 28.

(34) 3.5.3. Moving the Com mputation nal Mesh. Due to fluid flow the fibres of the oil filter f materrial are dissplaced. He ence for ph hysically ble results the t topolog gy of the fluid f mesh is adjusted d to the m moving solid d object plausib Recalcu ulation at each time e step is im mportant to t model the flow ffield chang ges and conservve geomettric consisttency. Both h regions, fluid and solid, s are h having a separate s mesh. For handlin ng both re egions a co ombined La agrangian and Eulerian grid me ethod is f region n a Lagran ngian grid is i employe ed. This me esh is adju usted to utilized. For the fluid oving boun ndary whicch is updatted every time step of the transient sim mulation. the mo Furtherrmore the e conserva ation valu ues, like mass, mo omentum and enerrgy are transpo orted with the movin ng grid. Fo or the solid d region a Eulerian m mesh is use ed. This means that the grid g is fixed d in time and a space. The displacement iinformation n is still commu unicated to o the fluid mesh lead ding to a deformatio d on of the fluid mesh.. As the solid mesh m is nott moved, for f post prrocessing itt is adjuste ed to the ffluid mesh h by the “warp” function of o ParaView w [16].. Figure 3.8: Fluid d (blue) an nd solid (o orange) re egions possses their own indep pendent mesh The im mportant co onstraint iss that at the t interface the verrtices and faces hav ve to be equal and a match congruentt with the opposite o re egion.. 29.

(35) Figure 3.9: Numb ber of faces and verttices on bo oth sides off the interfface has to o be the same and a congru uent The moving messh stays topological t ly the sam me, no ne ew cells w will be created or existing g ones desstroyed. In n OpenFOA AM® there are two types t of m mesh manip pulation approaches available. For the case of o no topo ology chan nges and m motion of interval mesh points p it is called “dyn namicMotio onSolverFv vMesh” [29 9,30]. Further on in order o to solve th he mesh motion m equation for calculating the mesh points an a appropriate e solver is necessary. The two available e solvers are: -. d displaceme entLaplaciaan: for th his solver the final displacement of the e mesh c componen ts is neede ed as well as the me esh displaccement of the intern nal field. I is based It d on the La aplacian difffusivity an nd the cell displaceme ent.. -. v velocityLap placian: this solver deals d with the t boundary velocitties instead d of the f final motio ons. It is used if the magnitu ude of the e maximum m displacement is s small compared to th he size of the overall domain. It is based d on the La aplacian d diffusivity and the ce ell motion velocity. v. For. a. continuo ous. mesh h. movem ment. and d. small. displacement. value es. the. velocityyLaplacian solver is the t most convenient c and hence e used forr all simula ations in this wo ork.. 30.

(36) The mesh motion consists of four steps: •. The moving interface is detected as declared by the user.. •. Points and faces are detected on the moving interface.. •. Displacement of solid is mapped to the fluid region. •. the fluid mesh is moved. The distortion of the mesh results in topology changes of each single cell volume. For simulation accuracy a good quality mesh is vital. Hence to keep the topology as consistent as possible over time, a propagation of the deformation values of single cells to the overall cell-collective is carried out. This will ensure that not only a single cell has to take all the deformation. It is dispersed in the overall collective. The movement of the boundary is predetermined by the boundary conditions calculated from fluid flow. The movement inside the fluid mesh and hence the coordinates of its grid points is controlled by a diffusion mechanism. It is based on the Laplace equation:. ∇ ⋅ (D f ⋅ ∇ u ) = 0. 3.27. where u represents the grid propagation velocity and Df the diffusion coefficient. Therefore the propagation of the deformation is ruled by the diffusion coefficient Df. The simplest way of improving mesh quality is by introducing variable diffusivity. There is a variety of coefficients available in OpenFOAM® [10,29,30] Distance-based methods: A number of boundary patches are selected by the user. The diffusion field Df is a function of cell centre distance L to the nearest selected boundary. -. Linear inverse Distance: Df = 1/L. -. Quadratic inverseDistance: Df= 1/L². -. Exponential: Df = 1/e(L). Quality-based methods: Here the diffusion field Df is a function of a cell quality measure: -. Mean cell non-orthogonality. -. Mean cell skewness 31.

(37) -. Mixed, wh hich is a co ombination n of mean cell non-o orthogonality and me ean cell s skewness.. All of the available diffusion n coefficien nts were te ested on a 2D mesh h with an arbitrary a square moving th hrough the computational mesh h.. Figure 3.10: A square moving from le eft to right through a 2D grid In Figure 3.10 on n the left hand side th he starting g point of the moving g square is shown. On the right hand d side it iss already in n motion. It is move ed to the m maximum possible p displace ement, wh here the siimulation is aborted due to he eavily disto orted grid cells. c At the tip of the squ uare, an area of inte erest is developing. In I Figure 3 3.11. it is marked with a red circle.. Figure 3.11: Both methodss not consserve the overall ge eometric consistency y of the nd side: “direction nal” diffussion coeffficient, riight hand d side: mesh (left han “faceOrrthogonalitty” diffusio on coefficie ent) On the e left hand side a quality ba ased diffussion coeffficient is u used. It iss called “directional” and takes the face skew wness into account. It I can be o observed that t the 32.

(38) cells att the tip arre heavily distorted and a even create c neg gative cell vvolumes. This T will quickly lead to te ermination of the sim mulation. The T diffusiion coefficcient, used on the and side, takes t the face f orthog gonality intto accountt. It showss that in im mminent right ha vicinity of the tip,, the cell volumes v arre staying in i good co ondition, whereas the e grid in p of the e fluid volu ume is disto orted. Again, this difffusion mecchanism is not the other parts best choice. Aftter testing g all diffussion coeffficients pro ovided byy OpenFOA AM®, a particular one sho owed exce ellent mesh h conservattion over time.. Figure 3.12: Favo ourable difffusion mecchanism: “q quadratic inverseDist i tance” Contrarry to the coefficientts used in Figure 3..11, in Fig gure 3.12 a distance e-based method d is used. It is called “quadrratic inversseDistance e” coefficie ent regard ding the fluid-so olid interfa ace. This diffusion mechanism m m was use ed for all of the fo ollowing simulattions. i is vital to t investiga ate the un nderlying diffusion d From this investigation follows that it mechan nism to gu uarantee the t conserrvation of the overa all mesh q quality durring the whole simulation s process.. 3.6. Structu ure of an n OpenF FOAM® Case C forr FSI Sim mulation n. In Ope enFOAM® the t data in n and outp put is hand dled over te ext files, ccalled dictio onaries. This ch hapter give es an overrview of th he necessa ary dictionaries and input para ameters for a su uccessful FSI F simulattion. Everyy case conssists of thre ee main diirectories. In each of them m there are e two sub folders, on ne for the fluid f and one o for the e solid region. The overall structure of o a typica al case for FSI simulation is sho own in Figu ure 3.13.. 33.

(39) Figure 3.13: Structure of a regularr OpenFOA AM® case with two o regions for FSI simulattion 1. “0”: This folder f conttains all the informattion about the initial conditions at time s step 0. F the fluid region (region1): For ( 1.) p: fo or initial co onditions of o the presssure (see chapter c 3.7 7.1) 2.) U: for f initial co onditions of o the veloccity (see chapter 3.7 7.1) 3.) cellM lMotionU/po pointMotion nU: for me esh movement (neccessary fo or mesh mottion solver: velocityLa aplacian) 4.) cellD lDisplaceme ment/pointD Displaceme ent: for mesh m move ement (for mesh mottion solver: displacem mentLaplaccian); not used u in thiss work F the sollid region (region2): For ( 1.) D: initial boun ndary conditions for the t displacement (see e chapter 3.7.1 3 ) 34.

(40) 2. “constant”: This directory comprises information about material properties, the computational mesh and moving mesh parameters. For the fluid region (region1): 1.) polyMesh: This folder contains a full description of the fluid mesh, like point coordinates, cells and boundaries. The text files are not supposed to be changed, as it would lead to severe inconsistencies and failures during simulation; 2.) transportProperties: This dictionary is used for adjusting fluid properties like the fluid density and the kinematic viscosity. 3.) fluidStructCouplingDict: This newly developed dictionary handles additional parameters like the ones for the collision concept. It gives the opportunity to define the fraction of the distance between fibres at which the collision concept shall start to work and to adjust the minimal distance for freezing of fibres. These topics will be discussed in detail in chapter 3.8. Further on it contains a switch to either turn the structure deformation on or off in case only a fluid flow solution is wanted by the user. This is also valid for the collision concept. 4.) dynamicMeshDict: This dictionary is responsible for all parameters concerning the mesh movement. It gives the user the opportunity to choose the diffusion mechanism, the necessary library and solver used for mesh movement. The necessary library for mesh movement without topology changes is called “libDynamicFvMesh.so” and is used for all simulations done in this work. For the solid region (region2): 1.) polyMesh: This folder contains a full description of the solid mesh. Again, the text files for points, cells and boundaries shall not be altered. 2.) mechanicalProperties: In this text file it is possible to adjust the properties of the solid, like the Young’s modulus, the Poisson ratio and the density of the material and if necessary neglect or consider plane stress. 35.

(41) 3.) couplingParameters: The user can define which boundary patch will be. the coupling interface between the solid and the fluid region. It must have the same name defined in each region in the boundary file of the polyMesh folder, which is predefined in the pre-processing of the mesh. With the parameter meshMotionTime the start time for mesh motion can be defined. Before this time point, only the fluid flow is calculated for stabilization of the fluid solution. The parameter motionRelaxation-iTime is used for relaxation of the solid displacement to support the convergence of the solution. 3. “system”: In this folder all the parameters for the solution procedure are handled. 1.) controlDict: Here all control parameters are set, for example time step size, start/end time of the overall simulation and the write accuracy. 2.) fvSchemes and fvSolution: Those dictionaries are used to set discretisation schemes, equation solvers, tolerances and all other algorithm controls for the run separately for each domain, region1 and region2. 3.) decomposeParDict:. In. this. dictionary,. the. method. for. domain. decomposition and its necessary coefficients for running the simulation in parallel are defined. More details on this topic will be discussed in chapter 4.1.1. The initial boundary field dictionaries for p and U are shown in chapter 3.7. A detailed description and further specifications of the system folder, with all available settings for discretisation schemes and solution controls are listed in the OpenFOAM® user’s and programmer’s guide [29,30]. The following figures show a selection of input dictionaries important for FSI simulation.. 36.

(42) Figure 3.14: Dictiionaries for mesh mo ovement. Figure 3.15: Dicttionaries fo or handling g domain properties: p transporttProperties for the omain and mechanica alPropertie es for the solid s domain fluid do. eshDict) a Figure 3.16: Dictionaries fo or mesh movement (dynamicM ( and auxiliarry input eters (fluidStructCoup plingDict) parame. 37.

(43) Figure 3.17: Dictiionary for general g co oupling parameters. 3.7. Applica ation to o a Simp ple Filterr Model. All com mponents of o the algo orithm werre assemb bled. The result r was a transien nt fluidstructure interacction solve er for inccompressib ble, lamin nar flow a and elastiic solid displace ement. In order to test t the de erived and programm med algoritthm a simp ple case was se et up. Forr the sake e of simpliicity this test t case consisted of a sing gle fibre situated d in a flow w channel at a right ang gle to the fluid f flow direction. d Materia al paramete ers and dim mensions of o the mod del:. Fluid:. Fiber:. Visccosity:. 1.76 1 x 10-4 m²/s. Den nsity:. 800 8 kg/m³. Leng gth of flow w channel L: L. 0.7 0 m. Crosss section A:. 0.16 0 m². Young’s modu ulus:. 10 1 5 N/m². gth of the fibre: Leng. 0.4 0 m. Diam meter of th he cross se ection:. 0.07 0 m. Figure 3.18: Dime ensioned fluid f domaiin with a siingle fibre at initial sttate 38.

(44) 3.7.1. Boundary conditions. There are two base types of boundaries used: 1. patch: this is the basic type, which contains no geometric or topology information about the mesh and hence is a generic patch. Here it is used for inlet and outlet. 2. wall: this is a type of patch including all the information about geometry and. topology of the defined patch. It is used for the surface of the fibres and the outer boundaries of the fluid region, besides inlet and outlet. The conditions applied on the patches are [30]: 1.) fixedValue: a Dirichlet boundary condition, where a certain value φ is specified; 2.) fixedGradient: a Neumann boundary condition, where the normal gradient of the value φ is specified; 3.) zeroGradient: a special type of the fixedGradient condition; the normal gradient of the value φ is set to zero; 4.) slip: for simulation of frictionless flow parallel to the patch; if the value φ is a scalar, it is set to zeroGradient, if the value φ is a vector the normal component and the normal gradient of the tangential component are set to zero. For this case the boundary conditions from the “0” dictionary reads as follows:. 39.

(45) Figure 3.19: Boun ndary field for the flu uid region. Figure 3.20: Boun ndary field for the so olid region In the solid regio on, for the e patch “Fibre” a spe ecial bound dary condiition is use ed. It is called “tractionDi “ isplacemen nt” and ha andles the applicatio on of pressure and traction forces on o the solid side (see e Figure 3..20).. 3.7.2. Results. The ressults of ap pplying the newly devveloped so olver to a single s fibre e are prese ented in this secction. Figu ure 3.21 sh hows the deformed d fibre. f Red marked re egions are regions with hig gh deformation value es. Blue marked area as are regions with zzero displaccement,. 40.

(46) as in th his case arre the end d points off the fibre, which are e kept at ffixed positions. In betwee en the disp placement is i gradual.. under fluid Figure 3.21: A sin ngle fibre deforming d d flow This firrst applica ation show ws very pla ausible ressults. In order o to vverify it, a simple analyticcal model was set up. In Figure 3.22 the deflectio on (displaccement) alo ong the length of the fibre is plotte ed for the e results acchieved byy simulatio on with the e newly ped FSI so olver and th he analyticcal model. develop. Figure 3.22: Co omparison of the deflection d (displacem ment) resu ults achieved by simulattion and th he ones callculated byy the analy ytical mode el It can be observe ed that the e two curvves match very well.. Especiallyy in the middle of the fibrre. The ma aximum diffference ap ppears closse to the fibre endpoints.. 41.

(47) For further testin ng, the casse was exte ended to six s fibres with w two different dia ameters (0.04 and a 0.06 m) m of the crross sectio on. The sam me model parameterrs and dim mensions as in th he previouss case werre used.. Figure 3.23: Model of six fib bres with different d th hickness be ending in th he flow of fluid As it ca an be seen n in Figure 3.23 the fibres are deformed due to the e flow of the t fluid in depe endency on n their thicckness. Fib bres with larger diam meters are barely defformed. In the middle of o the flow w channell, the flow w velocitie es of the fluid are at the maximu um. There efore fibre es with sm mall diame eters and additionally placed in the middle of the flow w channel are deform med the most. m Hence e not only the diame eter, but also the e spatial lo ocation of the t fibre in nfluences displaceme d nt. The flo ow of the flluid is deviiated by th he presence of the fib bres, which h result in regions of low or even ze ero velocitties. Corressponding to t this the absolute fluid flow velocity around the fibress is displaye ed in Figurre 3.24.. Figure 3.24: Velo ocity magniitude within surround ding fluid in n a cut pla ane 42.

(48) In this figure the e regions with w almosst zero ve elocity are marked blue. As exxpected, they ca an be found in the do ownstream m regions of o the fibre es. This is also confirrmed by the velocity vecttor field, shown s in the imme ediate vicin nity of the e fibres, which w is depicte ed in Figure e 3.25.. Figure 3.25: Velo ocity field of o fluid is deviated aro ound fibress Therefo ore, thin fiibres in the e downstre eam flow field f of thiicker fibress will also tend to zero de eformation, as there is no force e from the fluid flow acting a on tthe fibres. The forrces transfferred from m the fluid to the soliid are presssure and ttraction forces, as discusssed in chap pter 3.4. To o verify their implem mentation in nto the sollver, these e can be displayed. In Fig gure 3.26 on the left hand side the traction fforce varia ation is displayed, wherea as on the right r hand side the pressure p disstribution iis illustrate ed.. Figure 3.26: Tracction and pressure p ind duced by the t fluid. 43.

(49) It is evvident that pressure and tractio on forces mapped m on n the fibre es show ph hysically plausib ble results. The pressure is hig gh on the front side e of the fibres wherreas the traction n forces sh how high values v on the t sides. Therefore T el impleme entation their mode represe ents the exxpected ten ndencies co orrectly. The mo odel can no ow be exte ended to any a numbe er of artificiial fibres. F Figure 3.27 7 shows a case with sevveral fibress, randomly distributed and of o varying diameterss being med due to fluid flow. deform. Figure 3.27: Rand domly distrributed fibres in a flo ow channell. 3.8. Furthe er Impro ovement. 3.8.1. New Boundary Condition C ns. Until no ow the en nds of the fibres in the t simula ations have e been fixe ed. This does not reflect reality. Siince just a microsco opic piece e of the fiilter is mo odelled, fib bres do continu ue outside e of the modelling domain. Hence th he fibres should ha ave the possibility of glidiing through the dom main. As the ere are tw wo regions, one for so olid and one forr fluid, for both a new w boundaryy condition n has to be e applied. At first the bound dary condittion for the e flow velo ocity at the e walls of th he flow channel is change ed from the e fixed value of zero o to slip co ondition. Th he slip con ndition is provided p by Ope enFOAM® and a therefo ore it only needs to be b changed d inside the e “0” folde er of the case fo or the boun ndary cond dition of the e walls for the fluid.. 44.

(50) Figure 3.28: Dictiionary for the t fluid re egion with the new slip bounda ary conditio on To com mpare thosse two, th he fluid flo ow velocity y distributtion of the e cross se ection is displayed. In thiss case it is cut throug gh the firsst fiber in the t flow ch hannel, as can be seen in n Figure 3.2 29.. Figure 3.29: Fluid d flow velo ocity distribution in the t cross section s of the flow channel. c Bounda ary conditio on on the wall: w Left hand h side: fixedValue e 0, right h hand side: slip The diffference off those bou undary con nditions can clearly be b observe ed by zoom ming out the are ea of interrest, as dissplayed in Figure 3..30. The gap g shows the place e of the fibre, which w is cro ossing the area.. 45.

(51) Figure 3.30: Change of bou undary con ndition at the t wall fo or fluid flow w from fixe edValue = 0 to slip condittion On the e left hand side the previous p used bound dary condittion is sho own. It refe ers to a fixed va alue of zerro for the flow veloccity on the outer walll of the ch hannel, ma arked in blue. On O the righ ht hand sid de the slip condition is presente ed. “Slip” m means fricctionless flow, which w implie es that the e normal component of the vellocity vecto or and the normal gradien nt of the ta angential component are set to zero. Furtherr on, in ord der to mim mic reality in n a better way a new w boundaryy condition n for the solid re egion was implementted, in detail for the end pointss of the fib bres, which h are at contactt to the wa alls of the flow f chann nel. The ne ew conditio on means tthat the fib bres are not fixe ed anymorre but are e now able e to slide on their ends, e as th he fibres continue c over the walls of the siimulated domain. d Again, A thiss boundarry conditio on was develop ped especially for th his thesis as a it is not provided d by OpenFOAM®. The T idea behind this new boundary b c condition iss sketched in Figure 3.31.. Figure 3.31: Skettch of the new n bound dary condittion for the e end pointts of the fibres 46.

(52) Each fibre can slide along the t bound daries of th he computa ational dom main until a force retains it. This sp pring force Fs, which acts a reverssely on its motion, is defined ass. Fs = − k ⋅ Δx. 3.28. where k is the spring s constant, which is set to an emp pirical valu ue, and ∆xx is the sliding length disstance of th he fibre along the bo oundary. The T diction nary entry for this new bo oundary co ondition forr the solid region is shown s in th he following figure.. Figure 3.32: Dictiionary entrry for the solid s region n with new w boundaryy condition At the beginning of the mo ovement, the t fibre can slide un nhindered along the wall of w channel. This implies that th he force in ncreases with w distancce and eve entually the flow hinderss the fibre to move fu urther. Afte er some time, equilib brium deve elops betw ween the fluid fo orces acting g on the fiibre and th he spring force. f Whe en this point is reach hed, the fibre en nd points are regard ded as fixe ed and hence canno ot move fu urther, whe ereas in the mid ddle of the e fibre it ca an keep on moving an nd therefo ore will ben nd under th he force due to fluid flow. This beha aviour is sh hown in the e following figure.. 47.

(53) Figure 3.33: Sing gle fibre with fixed bo oundary co ondition off the fibre end pointss on the left han nd side and d spring bo oundary co ondition on n the right side. A clearr difference e can also o be observved in the e next Figu ure. The sa ame geom metry as displayed in Figurre 3.23 wa as used.. Figure 3.34: Diffference be etween fixe ed bounda ary conditio on in the upper figu ure and spring boundary condition c in the lowe er figure. 48.