## Publications

**Journal papers**

**Interface spaces for the Multiscale Robin Coupled Method in reservoir simulation**

GUIRALDELLO, R. T.; AUSAS, R. F.; **SOUSA, F. S.**; PEREIRA, F.; BUSCAGLIA, G. C.*Mathematics and Computers in Simulation*, 2018.

**Abstract:** The Multiscale Robin Coupled Method (MRCM) is a recent multiscale numerical method based on a non-overlapping domain decomposition procedure. One of its hallmarks is that the MRCM allows for the independent definition of interface spaces for pressure and flux over the skeleton of the decomposition. The accuracy of the MRCM depends on the choice of these interface spaces, as well as on an algorithmic parameter β in the Robin interface conditions imposed at the subdomain boundaries. This work presents an extensive numerical assessment of the MRCM in both of these aspects. Two types of interface spaces are implemented: usual piecewise polynomial spaces and informed spaces, the latter obtained from sets of snapshots by dimensionality reduction. Different distributions of the unknowns between pressure and flux are explored. Two nondimensionalizations of β are tested. The assessment is conducted on realistic, high contrast, channelized permeability fields from a SPE benchmark database. The results show that β, suitably nondimensionalized, can be fixed to unity to avoid any indeterminacy in the method. Further, with both types of spaces it is observed that a balanced distribution of the interface unknowns between pressure and flux renders the MRCM quite attractive both in accuracy and in computational cost, competitive with other multiscale methods from the literature.

**The Multiscale Robin Coupled Method for flows in porous media**

GUIRALDELLO, R. T.; AUSAS, R. F.; **SOUSA, F. S.**; PEREIRA, F.; BUSCAGLIA, G. C.*Journal of Computational Physics*, v. 355, p. 1-21, 2018.

**Abstract:** A multiscale mixed method aiming at the accurate approximation of velocity and pressure fields in heterogeneous porous media is proposed. The procedure is based on a new domain decomposition method in which the local problems are subject to Robin boundary conditions. The domain decomposition procedure is defined in terms of two independent spaces on the skeleton of the decomposition, corresponding to interface pressures and fluxes, that can be chosen with great flexibility to accommodate local features of the underlying permeability fields. The well-posedness of the new domain decomposition procedure is established and its connection with the method of Douglas et al. (1993) [12], is identified, also allowing us to reinterpret the known procedure as an optimized Schwarz (or Two-Lagrange-Multiplier) method. The multiscale property of the new domain decomposition method is indicated, and its relation with the Multiscale Mortar Mixed Finite Element Method (MMMFEM) and the Multiscale Hybrid-Mixed (MHM) Finite Element Method is discussed. Numerical simulations are presented aiming at illustrating several features of the new method. Initially we illustrate the possibility of switching from MMMFEM to MHM by suitably varying the Robin condition parameter in the new multiscale method. Then we turn our attention to realistic flows in high-contrast, channelized porous formations. We show that for a range of values of the Robin condition parameter our method provides better approximations for pressure and velocity than those computed with either the MMMFEM and the MHM. This is an indication that our method has the potential to produce more accurate velocity fields in the presence of rough, realistic permeability fields of petroleum reservoirs.

**Temporal large-eddy simulations of the lid-driven cavity by finite volume method**

CORREA, L.; MOMPEAN, G.; KUROKAWA, F. A.; **SOUSA, F. S.***Journal of the Brazilian Society of Mechanical Sciences and Engineering*, v. 40, p. 417, 2018.

**Abstract:** This paper describes in detail a numerical scheme to predict complex turbulent flows using a recent model based on temporal large-eddy simulations (TLES). To solve the equations a second-order finite volume numerical method coupled with a second-order time integration scheme is used. The numerical scheme is validated and then applied to present new results concerning the prediction of the complex turbulent flow in a cubic lid-driven cavity, at Reynolds numbers 𝑅𝑒=12,000 and 𝑅𝑒=18,000. The results obtained with the TLES are compared with direct numerical simulations and experimental data for the mean velocity flow field and for the Reynolds stresses, showing to be very attractive when compared to large-eddy simulations.

**Arbitrary Lagrangian-Eulerian (ALE)-based Finite Element Methods for rigid solids immersed in fluids**

PAZ, S.; **SOUSA, F. S.**; BUSCAGLIA, G. C. *Mecánica Computacional*, v. 34, p. 2165-2184, 2016.

**Abstract:** Arbitrary Lagrangian-Eulerian approaches are widely used in CFD, especially in multi-physics problems. They involve two tasks, namely the computation of the physical variables (velocity, stress, force, torque, etc.) and the determination of a suitable mesh deformation. We consider here a decoupled treatment of these two tasks, with high-order temporal schemes obtained by extrapolation, as discussed in F. Montefuscolo et.al. (J Comp Phys, 278:133-147, 2014) for capillary problems. Extensions of these schemes to fluid/rigid-body interaction are presented, adopting a variational formulation made popular by R. Glowinski et.al. in their work on Fictitious Domain Methods (J Comp Phys, 169:363- 426, 2001). The Arbitrary Lagrangian-Eulerian discretization turns the variational fluid-solid problem into a Differential-Algebraic Equation system for which several schemes, with different orders of accuracy, are implemented and evaluated. Special attention is dedicated to issues of stability, which is a fundamental obstacle towards the effective simulation of microfluidic fluid-solid interaction problems.

**Spurious transients of projection methods in micro flow simulations** **SOUSA, F. S.**; OISHI, C. M.; BUSCAGLIA, G. C. *Computer Methods in Applied Mechanics and Engineering*, v. 285, p. 659-693, 2015.

**Abstract:** The temporal behavior of projection methods for viscous incompressible low-Reynolds-number flows is addressed. The methods considered result from algebraically splitting the linear system corresponding to each time step, in such a way that the computation of velocity is segregated from that of pressure. Each method is characterized by two (possibly equal) approximate inverses (B_{1} and B_{2}) of the momentum-equation velocity matrix, plus a parameter γ which renders the method non-incremental (if γ=0) or incremental (if γ=1). The classical first-order projection method, together with more sophisticated methods (Perot’s second-order method, Yosida method, pseudo-exact factorization method) and their incremental variants are put into the same algebraic form and their accuracy numerically tested. Splitting errors of first, second and third order in the time step size δt are obtained, depending on the method. The methods are then discussed in terms of their ability and efficiency to compute steady states. Non-incremental methods are impractical because extremely small time steps are required for the steady state, which depends on δt, to be reasonably accurate. Incremental methods, on the other hand, either become unstable as δt is increased or develop a remarkable spurious transient which may last an extremely long time (much longer than any physical time scale involved). These transients have serious practical consequences on the simulation of steady (or slowly varying), low-inertia flows. From the physical viewpoint, the spurious transients may interfere with true slow processes of the system, such as heat transfer or species transport, without showing any obvious symptoms (wiggly behavior in space or time, for example, do not occur). From the computational viewpoint, the limitation in time step imposed by the spurious transient phenomenon weighs against choosing projection schemes for microflow applications, despite the low cost of each time step.

**High-order ALE schemes for incompressible capillary flows**

MONTEFUSCOLO, F.; **SOUSA, F. S.**; BUSCAGLIA, G. C. *Journal of Computational Physics*, v. 278, p. 133-147, 2014.

**Abstract:** The spatial discretization of problems with moving boundaries is considered, incorporating the temporal evolution of not just the mechanical variables, but also of the nodal positions of the moving mesh. The outcome is a system of Differential-Algebraic Equations (DAE) of index 2 or, in the case of inertialess flow, just 1. From the DAE formulation it its possible to define strategies to build schemes of arbitrary accuracy. We introduce here several schemes of order 2 and 3 that avoid the solution of a nonlinear system involving simultaneously the mechanical variables and the geometrical ones. This class of schemes has been the one adopted by the majority of practitioners of Computational Fluid Dynamics up to now. The proposed schemes indeed achieve the design accuracy, and further show stability restrictions that are not significantly more severe than those of popular first order schemes. The numerical experimentation is performed on capillary problems, discretized by both div-stable (P2/P1, P2/P1, P1+/P1) and equal-order (P1/P1, P1/P1 stabilized) finite elements, and incorporating surface tension and triple (contact) line effects.

**Uma técnica de correção de interface para o método 'Incompressible Smoothed Particle Hydrodynamics' **

CORDEIRO, D. F.; **SOUSA, F. S.**; CASTELO, A.; NÓBREGA, J. M. *Trends in Applied and Computational Mathematics*, v. 14(3), p. 347-358, 2013.

**Abstract:** The maintainability of a good particle distribution is of prior importance to SPH (Smoothed Particle Hydrodynamics) simulations, since it affects both stability and accuracy of the numerical results. When using ISPH (Incompressible SPH), the particle distribution can be improved by dislocating particles away from their Lagrangian trajectories. However, when applying this technique to multiphase flows, the correct position of the interface between fluids is disturbed, resulting in unphysical flow behavior. In this work, we propose a correction, based upon a constructed smooth distance function. Numerical tests show that this correction can recover the expected position of the interface. Finally, issues regarding mass conservation are assessed and reported.

**Numerical assessment of stability of interface discontinuous finite element pressure spaces** **SOUSA, F. S.**; AUSAS, R. F.; BUSCAGLIA, G. C. *Computer Methods in Applied Mechanics and Engineering*, v. 245-246, p. 63-74, 2012.

**Abstract:** The stability of two recently developed pressure spaces has been assessed numerically: The space proposed by Ausas *et al* (Comp. Meth. Appl. Mech. Eng., Vol. 199, 1019-1031, 2010), which is capable of representing discontinuous pressures, and the space proposed by Coppola-Owen and Codina (Int. J. Numer. Meth. Fluids, Vol. 49, 1287-1304, 2005), which can represent discontinuities in pressure gradients. We assess the stability of these spaces by numerically computing the inf-sup constants of several meshes. The inf-sup constant results as the solution of a generalized eigenvalue problems. Both spaces are in this way confirmed to be stable in their original form. An application of the same numerical assessment tool to the stabilized equal-order P_{1}/P_{1} formulation is then reported. An interesting finding is that the stabilization coefficient can be safely set to zero in an arbitrary band of elements without compromising the formulation's stability. An analogous result is also reported for the mini-element P_{1}^{+}/P_{1} when the velocity bubbles are removed in an arbitrary band of elements.

**An improved finite element space for discontinuous pressures**

AUSAS, R. F.; **SOUSA, F. S.**; BUSCAGLIA, G. C. *Computer Methods in Applied Mechanics and Engineering*, v. 199, p. 1019-1031, 2010.

**Abstract:** We consider incompressible Stokes flow with an internal interface at which the pressure is discontinuous, as happens for example in problems involving surface tension. We assume that the mesh does not follow the interface, which makes classical interpolation spaces to yield suboptimal convergence rates (typically, the interpolation error in the L2(Ω)-norm is of order h^1/2). We propose a modification of the P1-conforming space that accommodates discontinuities at the interface without introducing additional degrees of freedom or modifying the sparsity pattern of the linear system. The unknowns are the pressure values at the vertices of the mesh and the basis functions are computed locally at each element, so that the implementation of the proposed space into existing codes is straightforward. With this modification, numerical tests show that the interpolation order improves to O(h3/2). The new pressure space is implemented for the stable P1+/P1 source mini-element discretization, and for the stabilized equal-order P1/P1 discretization. Assessment is carried out for Poiseuille flow with a forcing surface and for a static bubble. In all cases the proposed pressure space leads to improved convergence orders and to more accurate results than the standard P1 space. In addition, two NavierStokes simulations with moving interfaces (RayleighTaylor instability and merging bubbles) are reported to show that the proposed space is robust enough to carry out realistic simulations.

**The MAC method**

MCKEE, S.; TOMÉ, M. F.; FERREIRA, V. G.; CUMINATO, J. A.; CASTELO, A.; **SOUSA, F. S.**; MANGIAVACCHI, N. *Computers & Fluids*, v. 37, p. 907-930, 2008.

**Abstract:** In this article recent advances in the Marker and Cell (MAC) method will be reviewed. The MAC technique dates back to the early 1960s at the Los Alamos Laboratories and this article starts with a historical review, and then a brief discussion of related techniques. Improvements since the early days of MAC (and the Simplified MAC - SMAC) include automatic time-stepping, the use of the conjugate gradient method to solve the Poisson equation for the corrected velocity potential, greater efficiency through stripping out the virtual particles (markers) other than those near the free surface, and more accurate approximations of the free surface boundary conditions, the addition of bounded high accuracy upwinding for the convected terms (thereby being able to solve higher Reynolds number flows), and a (dynamics) flow visualization facility. More recently, effective techniques for surface and interfacial flows and, in particular, for accurately tracking the associated surface(s)/interface(s) including moving contact angles have been developed. This article will concentrate principally on a three-dimensional version of the SMAC method. It will eschew both code verification and model validation; instead it will emphasize the applications that the MAC method can solve, from multiphase flows to rheology.

**Local volume-conserving free surface smoothing** **SOUSA, F. S.**; CASTELO, A.; NONATO, L. G.; MANGIAVACCHI, N.; CUMINATO, J. A. *Communications in Numerical Methods in Engineering*, v. 23, p. 109-120, 2007.

**Abstract:** Removing high-frequency undulations in surfaces is a problem that appears in different fields, such as computer graphics and computational fluid mechanics. This problem is typically handled by surface smoothing techniques, such as Laplacian filters, that eliminate high-frequency undulations but degrade the volume encompassed by the surface. The need for conserving volume (or mass) rules out the use of such techniques in several application, as for example incompressible flows. In this work we present a smoothing technique that suppresses undulations while conserving local volumes, ensuring that the global mass is conserved. The effectiveness of the proposed technique is illustrated in synthetic datasets as well as in free surface flows simulation. Comparisons between our smoothing approach and the well-known Laplacian filter are also presented.

**A Lagrangian level-set approach for the simulation of incompressible two-fluid flows** **SOUSA, F. S.**; MANGIAVACCHI, N.*International Journal for Numerical Methods in Fluids*, v. 47, p. 1393-1401, 2005.

**Abstract:** A Lagrangian level-set method to solve incompressible two-dimensional two-fluid flows is presented. The Navier-Stokes equations are discretized by a Galerkin finite element method. A projection method is employed to decouple the system of non-linear equations. The interface between fluids is represented by the zero level set of a function plus additional marker points of the computational mesh. In the standard Eulerian level-set method, this function is advected through the domain by solving a pure advection equation. To reduce mass conservation errors that can arise from this step, our method employs a Lagrangian technique which moves the nodes of the finite element mesh, and consequently, the information stored in each node. The quality of the mesh is controlled by a remeshing procedure, avoiding bad triangles by flipping edges, inserting or removing vertices from the triangulation. Results of numerical simulations are presented, illustrating the improvements in mass conservation and accuracy of this new methodology.

**A numerical method for solving three-dimensional generalized Newtonian free surface flows**

TOMÉ, M. F.; GROSSI, L.; CASTELO, A.; CUMINATO, J. A.; MANGIAVACCHI, N.; FERREIRA, V. G.; **SOUSA, F. S.**; MCKEE, S. *Journal of Non-Newtonian Fluid Mechanics*, v. 123, p. 85-103, 2004.

**Abstract:** This work presents a numerical technique for solving three-dimensional generalized Newtonian free surface flows. It is an extension to three dimensions of the technique introduced by Tomé et al. [M.F. Tomé, B. Duffy, S. McKee, A numerical technique for solving unsteady non-Newtonian free surface flows, J. Non-Newtonian Fluid Mech. 62 (1996) 934] but additionally includes many other features. The governing equations are solved by a finite difference method on a staggered grid. It uses marker particles to describe the fluid; these particles provide the location and visualization of the fluid free surface. As currently implemented, the present method can simulate generalized Newtonian flow in which the viscosity is modelled using the Cross model. The numerical technique presented in this paper is validated by using exact solutions for the flow of a Cross model fluid inside a pipe and convergence is demonstrated by means of grid refinement for the problem of a spreading drop. Numerical results showing the flow of a generalized Newtonian fluid jet impinging onto a flat surface and that of a jet buckling are given.

**A front-tracking/front-capturing method for the simulation of 3D multi-fluid flows with free surfaces** **SOUSA, F. S.**; MANGIAVACCHI, N.; NONATO, L. G.; CASTELO, A.; TOMÉ, M. F.; FERREIRA, V. G.; CUMINATO, J. A.; MCKEE, S.*Journal of Computational Physics*, v. 198, p. 469-499, 2004.

**Abstract:** A method for simulating incompressible, imiscible, unsteady, Newtonian, multi-fluid flows with free surfaces is described. A sharp interface separates fluids of different density and viscosity. Surface and interfacial tensions are also considered and the required curvature is geometrically approximated at the fronts by a least squares quadratic fitting. To remove small undulations at the fronts, a mass-conserving filter is employed. The numerical method employed to solve the Navier-Stokes equations is based on the GENSMAC-3D front-tracking method. The velocity field is computed using a finite-difference scheme on an Eulerian grid. The free-surface and the interfaces are represented by an unstructured Lagrangian grid moving through an Eulerian grid. The method was validated by comparing the numerical results with analytical results for a number of simple problems. Complex numerical simulations show the capability and emphasize the robustness of this new method.