Acta Polytechnica doi:10.14311/AP.2017.57.0245 Acta Polytechnica 57(4):245–251, 2017 © Czech Technical University in Prague, 2017 available online at http://ojs.cvut.cz/ojs/index.php/ap FICTITIOUS DOMAIN METHOD FOR NUMERICAL SIMULATION OF INCOMPRESSIBLE VISCOUS FLOW AROUND RIGID BODIES Matej Beňo∗, Bořek Patzák Czech Technical University in Prague, Faculty of Civil Engineering, Department of Mechanics, Thákurova 7, 166 29 Prague, Czech Republic ∗ corresponding author: matej.beno@fsv.cvut.cz Abstract. This article describes the method of efficient simulation of the flow around potentially many rigid obstacles. The finite element implementation is based on the incompressible Navier-Stokes equations using structured, regular, two dimensional triangular mesh. The fictitious domain method is introduced to account for the presence of rigid particles, representing obstacles to the flow. To enforce rigid body constraints in parts corresponding to rigid obstacles, Lagrange multipliers are used. For time discretization, an operator splitting technique is used. The model is validated using 2D channel flow simulations with circular obstacles. Different possibilities of enforcing rigid body constraints are compared to the fully resolved simulations and optimal strategy is recommended. Keywords: finite element method; computational fluid dynamics; fictitious domain method; Lagrange multipliers; flow around obstacles. 1. Introduction The motivation of the paper comes from the mod- eling of the fresh concrete casting. The simulation models of fresh concrete aiming at structural-scale ap- plications typically consider concrete suspension as a homogeneous, non-Newtonian fluid, whose rheological properties can be derived from the mix composition [1]. The flow can be then efficiently described by Navier- Stokes equations considering non-Newtonian fluid and solved using the Finite Element Method (FEM), for example. Of course, the efficiency of the homogeneous approach comes at the expense of the coarser descrip- tion of the flow, where sub-scale phenomena can only be approximately accounted for by post-processing simulation results, e.g. to determine the distribution and orientation of reinforcing fibers [2], or by heuris- tic modification of constitutive parameters, e.g. to account for the effect of traditional reinforcement [3]. Especially the latter aspect is critical in the modeling of casting processes in highly-reinforced structures that represent the major field of application for self- compacting concrete. The explicit representation of individual reinforcing bars in the FE model would lead to extremely fine computational meshes and will result in extremely high computational demands in terms of resources and time. This article aims at reducing these compu- tational costs by avoiding the need of explicit repre- sentation of individual reinforcing bars by adopting an approach based on fictitious domain method to solve the problem of Newtonian incompressible flow with rigid body obstacles on regular, structured computa- tional grid, where the individual reinforcing bars can be inserted arbitrarily and independently of the un- derlying mesh. The individual bars are accounted for by enforcing no-flow constraints in a volume occupied by the bar using Lagrange multipliers. To solve the incompressible Navier-Stokes problem, a 2D Eulerian formulation of finite element method is developed, using the “Taylor-Hood” P2/P1 elements, which are quadratic in velocity and linear in pressure, satisfy- ing the LBB condition [4]. The time discretization is based on operator splitting approach, in particular, the fractional-time-step scheme described by Marchuk has been employed [5]. The fundamental idea of fictitious domain method consists in extending the problem to one geometri- cally simpler domain covering both fluid and obstacles – the so-called fictitious domain. The no-flow con- straints in subdomains corresponding to individual bars are enforced using distributed Lagrange multipli- ers, which represent the additional body forces needed to enforce zero flow inside the regions representing obstacles. The pressure is constrained by incompress- ibility of the fluid in the fictitious domain. This so called distributed Lagrange multiplier fictitious do- main method has been introduced by Glowinski et al. [6]. The method has then been used by Bertrand et al. [7] and Tanguy et al. [8] to calculate three- dimensional flows. Recently, the method has been extended by Baaijens [9] and Yu [10] to handle fluid- structure interactions. The main advantage of this method is a much sim- pler generation of the computational grid, because the mesh does not need to conform to the geome- try of the reinforcement. It can be fully structured and regular, allowing to take advantage of specialized, fast numerical solvers. This feature can significantly reduce computational requirements of the problem. To the authors knowledge, there has been no study published investigating the influence of the sampling point selection on the quality of the solution. The 245 http://dx.doi.org/10.14311/AP.2017.57.0245 http://ojs.cvut.cz/ojs/index.php/ap Matej Beňo, Bořek Patzák Acta Polytechnica Figure 1. Original and fictitious domain problems. goal of this paper is to propose and evaluate differ- ent strategies for enforcing no-flow constraints using different sets of sampling points and compare their performance to the fully resolved simulations. 1.1. Fictitious domain method The basic idea of the method is based on a cancellation of the forces and moments between obstacles and fluid in the combined weak formulation for particle-fluid motion solved everywhere in the fictitious domain. To enforce the no-flow constraints (sometimes referred as rigid-body constraints when moving rigid particles are considered) to account for rigid particles, Lagrange multipliers are introduced. These multipliers represent the additional body forces needed to maintain the rigid body motion inside the moving particles or the no flow constraint inside the regions representing fixed particles. As already mentioned, the advantage of this approach consist in ability to use geometrically simpler, regular meshes. 2. Problem formulation 2.1. The governing equation of the flow Assume an incompressible, viscous, Newtonian fluid occupying at the given time t ∈ (0,T) the delimited domain Ω ⊂ R2 with boundary Γ. Let denote by x = {xi}2i=1 a generic point in Ω. Let further denote by u(x, t) velocity and by p(x, t) pressure, both governed by Navier-Stokes equations % du dt = %g + ∇·σ on Ω \P, (1) ∇·u = 0 on Ω \P, (2) where % is density of the fluid, u is velocity of the fluid, and σ is stress tensor. For an incompressible, Newtonian viscous fluid, the stress can be decomposed into hydrostatic and deviatoric components σ = −pI + 2ηD[u], (3) where p is hydrostatic pressure in the fluid, η is the viscosity (assumed constant), and 2ηD[u] is deviatoric stress tensor. Relations (1)–(3) are to be complemented by the appropriate initial and boundary conditions u(x, 0) = u0 on Ω \P, (4) ∇·u0 = 0, (5) u(x, t) = uΓ(t) on Γ, (6)∫ Γ uΓ(t)n̂dx = 0, (7) where n̂ is unit normal vector pointing out of Γ. 2.2. Governing equations for obstacles Hydrodynamic force Fi and torque Ti acting on the i-th particle can be evaluated by summing up corre- sponding fluid induced forces acting on the particle boundary: Fi = ∫ ∂Pi σn̂ ds, Ti = ∫ ∂Pi rixσn̂ ds, (8) where Pi is the region occupied by the particle and ri = x − Xi is relative position vector to particle center Xi. 3. Weak form Let us introduce the approximation of spaces for test and trial functions for velocities and pressure in the fluid part of domain Ω \P without obstacles: W = {u ∈ H1(Ω \P) : u = uΓ(t) on Γ}, V = {v ∈ H1(Ω \P) : v = 0 on Γ}, L = {p ∈ L2(Ω \P)}, K = {q ∈ L2(Ω \P)}. By using the method of weighted residuals applied to (1), (2) and using finite dimensional spaces W0.h, V0.h, L0.h, K0.h, approximating W, V, L, K defined above we arrive at the following finite-element approx- imation of the Navier-Stokes equations: Find uh ∈ W0.h, p ∈ L0.h satisfying∫ Ω\P ( % ∂uh ∂t + (uh ·∇)uh ) ·vh dx − ∫ Ω\P ph∇·vh dx + ∫ Ω\P 2ηD[uh] : D[vh] dx = 0 for all vh ∈ W0.h, (9)∫ Ω\P qh∇·uh dx = 0 for all qh ∈ L0.h, (10) uh(0) = u0,h on Ω, (11) where u0,h is a divergence-free initial velocity. Since in (9) u is divergence-free and satisfies the Dirichlet boundary condition (4) on Γ, we can write ∫ Ω\P 2ηD[uh] : D[vh] dx = ∫ Ω\P η∇uh : ∇vh dx for all vh ∈ W0.h. 246 vol. 57 no. 4/2017 Fictitious Domain Method for a Numerical Simulation 3.1. A fictitious domain formulation To extend from the domain Ω \P to the fictitious domain Ω, it is necessary to enforce no-flow constraints inside each Pi. Let us introduce an approximation of spaces for test and trial functions velocities and pressure for the fictitious domain Ω: W = {u ∈ H1(Ω) : u = uΓ(t) on Γ}, V = {v ∈ H1(Ω) : v = 0 on Γ}, L = {p ∈ L2(Ω)}, K = {q ∈ L2(Ω)}. By extending the incompressibility condition to the whole Ω, we get∫ Ω q∇·u dx = 0 for all q ∈ L2h(Ω). The condition to enforce rigid body (or no-flow) con- strains for each obstacle can be expressed as u(x, t) −u(Xi, t) = 0 for all x ∈ P. (12) To relax the no-flow constraint defined above, the family of Lagrange multipliers λj is introduced rep- resenting a discrete set of points covering Pj, such that Λh = { µh :µh = M∑ i=1 µiδ(x−xi), µ1,...,µm ∈R2 } , where δ is the Dirac delta function at x = 0. Using the space Λh defined above, the condition (12) is relaxed to〈 µh,u(x, t) −u(Xi, t) 〉 = 0 for all µh ∈ Λh, (13) where the inner product on the obstacle is defined as 〈µh,vh〉P = M∑ i=1 µh,i ·vh(xi). In case of a fixed particle, it is necessary to prescribe zero velocity at the particle center. In case of moving particle, the equations of motions for a particle have to be solved as well. The combined weak solution is then to find uh ∈ W0.h, ph ∈ L0.h, where finite dimensional spaces W0.h, L0.h are approximating W, L spaces defined above and λh ∈ Λh is such that∫ Ω % (∂u ∂t + (uh ·∇)uh ) vh dx− ∫ Ω ph∇·vh dx + ∫ Ω 2ηD[uh] : D[vh] dx = 〈λh,vh〉P for all vh ∈ Vh,∫ Ω qh∇·uh dx = 0 for all qh ∈ L2h, 〈µh,uh〉P = 0 for all µh ∈ Λh, u(0) = u0 on Ω. 3.2. Time discretization by operator splitting As pointed out by Glowinski [11], numerical solutions of the nonlinear partial differential Navier-Stokes equa- tions are not trivial, mainly due to the following rea- sons: (i) the above equations are nonlinear and (ii) the incompressibility condition is difficult to handle. The above equations represent a system of partial differential equations, coupled through the nonlinear term , the incompressibility condition ∇·u = 0, and sometimes through the boundary conditions. The use of time discretization by operator split- ting will partly overcome the above difficulties; in particular, decoupling of difficulties associated with the nonlinearity from those associated with the incom- pressibility condition. Following the work of Glowin- ski and Pironneau [12], assuming the following initial value problem dϕ dt + A(ϕ) = 0, ϕ(0) = ϕ0, where A is an operator (possibly nonlinear, and even multivalued) from a Hilbert space H into itself and where ϕ0 ∈ H. A number of splitting techniques have been proposed, for an overview, we refer to [5]. In this work, we have used the scheme proposed by Marchuk [5], which is described in the next section. 3.3. Marchuk’s fractional-step scheme Marchuk assumed a decomposition of the operator A into the following nontrivial decomposition A = A1 + A2 + A3 (14) (by nontrivial, we mean that operators A1, A2 and A3 are individually simpler than A). By assuming a time discretization with the time step ∆t, the updated value at the end of the time step can be computed using a three step integration procedure as follows: ϕn+1/3 −ϕn ∆t + A1(ϕn+1/3) = fn+11 , ϕn+2/3 −ϕn+1/3 ∆t + A2(ϕn+2/3) = fn+12 , ϕn+1 −ϕn+2/3 ∆t + A3(ϕn+1) = fn+13 . By applying an operator splitting to the problem (9)– (12), one obtains (with 0 ≤ α,β ≤ 1 and α + β = 1): Find un+1/3h Wh and p n+1/3 ∈ Lh such that % ∫ Ω u n+1/3 h −u n h ∆t ·vh dx− ∫ Ω p n+1/3 h ∇·vh dx = 0 for all vh ∈ Vh, (15)∫ Ω qh∇·u n+1/3 h dx = 0 for all qh ∈ L 2 h. (16) 247 Matej Beňo, Bořek Patzák Acta Polytechnica Figure 2. The different strategies for selection of sampling points. The boxed cross represents boundary condition application (zero velocity) in case of fixed particle, the open cross indicate location of sampling point. Find un+2/3h ∈ Wh such that % ∫ Ω u n+2/3 h −u n+1/3 h ∆t ·vh dx −% ∫ Ω (un+1/3h ·∇)u n+2/3 h ·vh dx + 2αη ∫ Ω D[un+2/3h ] : D[vh] dx = 0 for all vh ∈ Vh. (17) Finally find un+1h ∈ Wh, λ n+1 h ∈ Λh such that % ∫ Ω un+1h −u n+2/3 h ∆t ·vh dx + 2βη ∫ Ω D[un+1h ] : D[vh] dx = 〈λ n+1 h ,vh〉 for all vh ∈ Vh, (18) 〈µh,un+1h 〉P = 0 for all µh ∈ Λh. (19) This scheme is only first order accurate, but with good stability properties, see [5]. The described method has also been used in [6]. 4. Sampling point selection As already mentioned, the rigid body or no flow con- straints are imposed in a weak sense using Lagrange multipliers. In a practical implementation, the corre- sponding integral in a weak form is evaluated as a sum over a set of discrete sampling points. In this paper, we have proposed and tested several strategies to en- force constrains by combining one or more Lagrange multipliers per particle and using different strategies for the sampling point selection. For the considered Figure 3. The problem of 2D flow around single obstacle. Figure 4. The mesh of 2D flow with single obstacle. circular geometry of particle, the proposed strategies consist of generating the different patterns of sam- pling points in radial and longitudinal directions, see Figure 2. These results of these strategies are com- pared using two test scenarios to results obtained with fully resolved simulations, where the geometry of each individual particle is represented exactly by the mesh and corresponding no-flow constraints are exactly satisfied. 5. Numerical tests 5.1. Flow around single obstacle The first example to evaluate proposed strategies solves the 2D flow around a single, rigid particle. The geometry of the problem is depicted in Fig- ure 3, together with initial and boundary conditions. The perfect friction (no slip) boundary condition has been assumed on horizontal edges. The domain has been discretized into 1820 triangular Taylor-Hood ele- ments, see Figure 4 top. The values of mass density % = 1.0 kg/m3, viscosity η = 10−1 Pa s, and time step ∆t = 0.001 s have been used in this study. The described model has been implemented into 248 vol. 57 no. 4/2017 Fictitious Domain Method for a Numerical Simulation Figure 5. Velocity profiles through obstacle. Figure 6. The geometry and the boundary conditions of the flow in tube with obstacle (448 el.). an object oriented FE code oofem [13]. The differ- ent strategies for a selection of sampling points have been evaluated and compared to results obtained by fully resolved simulation (see Figure 4 bottom for discretization, consisting of 1656 elements and 3462 nodes). The comparison with fully resolved simulation has been made at a steady state. The individual strategies considered are described in Table 1. The Figure 5 shows the profiles of velocity (horizontal component) in a vertical section passing through the particle center. Table 1 also contains the L2 norm of the difference of the obtained velocity and velocity from fully resolved simulation along the vertical profile passing the center of the particle. The obtained results indicate that an optimum strategy for sampling point selection is to have multiple rings with a dedicated multiplier for each ring. The best result (in terms of smaller difference) has been obtained for strategy 05 consisting of 12 rings and 12 multipliers for each ring. However, it is necessary to balance the computational effort (related to number of rings and number of multipliers) and obtained error. From this point of view, the optimal seems to be the strategy 03 with 4 rings. 5.2. Flow around two obstacles The second example represents a flow in a tube with two particles, as illustrated in Figure 6. In this test, three different discretizations of the fictitious domain have been used with an increasing mesh density, con- taining 448, 1792, and 7168 elements. The perfect friction on both horizontal edges (no slip condition) has been assumed. The results are again compared by means of comparing the velocity profiles along the 249 Matej Beňo, Bořek Patzák Acta Polytechnica No. Lagrange multipliers Figure Error L2-norm 10 fully resolved simulation 0 01 1 multiplier per ring – 1 ring 15 A 58.10 02 1 multiplier per ring – 2 rings 15 B 9.36 03 1 multiplier per ring – 4 rings 15 B 1.13 04 1 multiplier per ring – 8 rings 15 B 0.93 05 1 multiplier per ring – 12 rings 15 B 0.89 001 1 multiplier for 2 rings 15 C 22.66 002 1 multiplier for 4 rings 15 C 7.48 003 1 multiplier for 8 rings 15 C 4.45 004 1 multiplier for 12 rings 15 C 4.04 01-P1 same as 01 – sporadic points 15 D 56.33 02-P1 same as 02 – sporadic points 15 D 21.10 03-P1 same as 03 – sporadic points 15 D 13.92 04-P1 same as 04 – sporadic points 15 D 9.76 05-P1 same as 05 – sporadic points 15 D 8.88 Table 1. Different ways of introducing of sample points with L2 norm of velocity difference in vertical profile passing the center of obstacle. Figure 7. Vertical profiles of velocity passing the centers of both particles. vertical section, indicated in Figure 6. Four (eight, twenty-four) Lagrange multipliers per obstacle hav been used. Figure 6 shows obtained velocity profiles using the fictitious domain method and using the fully resolved simulation. The obtained results show excellent agreement to the results obtained by the fully resolved simulation. 6. Conclusion The presented paper deals with the modelling of incom- pressible flow around fixed (or moving, rigid) particles using the concept of fictitious domain method and demonstrates the capability of this method to describe the flow with solid obstacles. The paper demonstrates that the accuracy of the method depends on the strategy of choosing sampling points in addition to traditional aspects, such as time stepping and integration or space discretization. From the simulations performed using different strategies, the uniform distribution of sampling points is strongly recommended, while each element should contain at least 1 sampling point. Although the obtained results are problem and grid dependent, they demonstrate the clear dependence of the sampling point selection on the quality of the solution. Acknowledgements The authors would like to acknowledge the support of the Czech Science Foundation under the project 13-23584S. 250 vol. 57 no. 4/2017 Fictitious Domain Method for a Numerical Simulation References [1] H. Mori, Y. Tanigawa. Simulation methods for fluidity of fresh concrete. Memoirs of the School of Engineering. 1992, Vol. 1, 44. [2] F. Kolařík, B. Patzák, L.N. Thrane. Modeling of fiber orientation in viscous fluid flow with application to self-compacting concrete. Computers & Structures. 2015, Vol. 154, pp. 91-100. [3] K. Vassilic, B. Meng, H.C. Kühne, N. Roussel. Flow of fresh concrete through steel bars: A porous medium analogy. Cement and Concrete Research. May 2011, Vol. 41, 5, pp. 496-503. [4] Babuška, I. The Finite Element Method with Lagrangian Multipliers. Numerische Mathematik. 20, pp. 179-192. [5] Marchuk, G.I. Splitting and Alternating Direction Methods. [book auth.] J.L. Lions P.G. Ciarlet. Handbook of Numerical Analysis. Amsterdam, North-Holland, 1990, Vol. 1, pp. 197-461. [6] R. Glowinski, T.-W. Pan, T.I. Hesla, D.D. Joseph. A distributed Lagrange multiplier/fictitious domain method for particulate flows. International Journal of Multiphase Flow. 1999, Vol. 25, pp. 755-794. [7] F. Bertrand, P.A. Tanguy, F. Thibault. A three-dimensional fictitious domain method for incompressible fluid flow problems. Int. J. Numer. Meth. Fluids. 1997, Vol. 25, pp. 719-736. [8] P.A. Tanguy, F. Bertrand, R. Labie, E. Brito- De La Fuente. Numerical modelling of the mixing of viscoplastic slurries in a twin-blade planetary mixer. Trans. IChemE. 1996, Vol. 74 (Part A), pp. 499-504. [9] Baaijens, F.P.T. A fictitious domain/mortar element method for fluid–structure interaction. Int. J. Numer. Meth. Fluids. 2001, Vol. 35, pp. 743-761. [10] Yu, Z. A DLM/FD method for fluid/flexible-body interactions. J. Comput. Phys. 2005, Vol. 207, pp. 1-27. [11] Glowinski, R. Finite Element Methods for Incompressible Viscous Flow. [book auth.] J.L. Lions P.G. Ciarlet. Handbook of Numerical Analysis. Amsterdam, North-Hollad, 2003, Vol. 9, pp. 3-76. [12] R. Glowinski, O. Pironneau. Finite element methods for Navier-Stokes equations. Annual Reviews Fluid Mech. 1992, Vol. 24, pp. 167-204. [13] Patzák, B. OOFEM- an object-oriented simulation tool for advanced modeling of materials and structures. Acta Polytechnica. 2012, Vol. 52, 6, pp. 59-66. 251 Acta Polytechnica 57(4):245–251, 2017 1 Introduction 1.1 Fictitious domain method 2 Problem formulation 2.1 The governing equation of the flow 2.2 Governing equations for obstacles 3 Weak form 3.1 A fictitious domain formulation 3.2 Time discretization by operator splitting 3.3 Marchuk’s fractional-step scheme 4 Sampling point selection 5 Numerical tests 5.1 Flow around single obstacle 5.2 Flow around two obstacles 6 Conclusion Acknowledgements References