An improved immersed boundary method by coupling of the multi-direct forcing and Fourier pseudo-spectral methods

The present paper aims to report an improve on the development of the IMERSPEC methodology, which is a methodology that couples the Fourier pseudo-spectral and immersed boundary methods, which can now solve flows over complex geometries. In the present work, the Lagrangian mesh that represents the immersed body inside the flow does not have to coincide with the Eulerian mesh. Furthermore, a high accuracy for the resultant methodology is achieved. These two important improvements arise from the development of the multi-direct forcing method and used to calculate the force field of the immersed boundary method. In order to verify the IMERSPEC methodology, a manufactured solution technique was used. The results of a refined mesh were employed, and different distribution and interpolation functions were used and analysed. Finally, flows over a backward-facing step geometry were simulated in two dimensions. This type of flow presents several physical features, for example, detachment and reattachment of the boundary layer and vortex generation, which are difficult to capture using numerical solutions. A comparison with the experimental data is presented and shows that the IMERSPEC methodology is promising for computational fluid dynamics applications and its analysis.


Introduction
Computational fluid dynamics (CFD) methodologies have became quite common and reliable for both engineering and scientific problems in academic and industrial environments. The developments achieved until this point have enabled the increase in their potential, as well as the level of confidence that the scientific community places in them. Large CFD applications have been conducted by several industrial sectors, such as petroleum, chemical and aeronautical industries. The common macrofeatures for industrial applications are physical and geometrical complexity, flexible and moving geometries, turbulent flows and multiple characteristic scales. In order to solve such problems, special methodologies have been used (i.e. adaptive meshes over complex geometries) to achieve high accuracy and a high rate of convergence. On the other hand, the computational cost increases for more accurate methodologies, both in CPU time and in storage memory. For example, in studies of aeroacoustics, combustion and surfactant transport, a low concentration is normally sufficient to model large physical effects. Thus, high-accuracy solutions for PDEs are necessary.
The current perception of researchers in the field of mathematical modelling is that the solution of differential models should have the following features: high accuracy, high convergence rates, low computational cost (memory storage and CPU time), ease in handling complex, moving and deformable geometries, ease in programming and parallelising and good scalability. The techniques used to solve flow problems over complex and moving geometries require moving and deforming meshes or even remaking the mesh after a given number of time steps. These procedures cause the methodologies to have a high computational cost and have a more complex implementation [1].
The immersed boundary methodology has been developed since 1970 by several researchers [2][3][4][5][6]. This method fulfils a portion of the requirements described above; specifically, it can handle complex and moving geometries. This methodology, in contrast to unstructured grids, allows flows to be simulated over complex geometries using only a single Cartesian mesh. On the other hand, finite volume, finite element and finite difference methodologies involve large computational costs to reach high-accuracy solutions [1]. Furthermore, using spectral methods like collocation, Galerkin, Tau, etc., high-accuracy solutions are achieved, see the works of [7][8][9][10].
Computer simulations with high accuracy can be performed by solving fluid flow equations with numerical methods of high order of convergence [9,11], e.g. a CFD problem is solved in a given mesh with a high-order convergence method (greater than two), and results are obtained with high accuracy, that is, with insignificant errors. At the same time, using a method of convergence order, less than two, requires a much finer mesh to achieve the same level of accuracy in the results of the same problem. Sometimes, this level of refinement makes it impossible to carry out a simulation, as it requires a greater amount of RAM and processing time.
Among the high-order numerical convergence methods, the spectral methods stand out, as they are characterised by an exponential convergence rate in relation to the mesh refinement, [7,10,[12][13][14]. Specifically, the collocation Fourier pseudo-spectral method is especially used in the solution of the Navier-Stokes equations, because it eliminates the pressure gradient solution for incompressible flows, by guaranteeing the conservation of mass through of the projection procedure of the nonlinear term. In mathematical and computational terms, the linear system solution (to solve the Poisson equation) is replaced by a matrix vector product, this is a numerical procedure faster to performed computationally and that demands less RAM memory consumption. In any case, the pressure field can be obtained in post-processing.
The main computational cost (processing time and RAM memory usage) when using the pseudo-spectral Fourier method is in the process of transforming the velocity fields into spectral space or Fourier space. For that, the algorithm called fast Fourier transform is used [15,16], which is still, sometimes, more efficient than the solution of linear systems, mainly, when comparing simulations with more mesh refinement. However, the pseudo-spectral Fourier collocation method requires periodic boundary conditions [17], which limits its use in computational fluid dynamics to few useful problems, [13,18].
The merging of immersed boundary and Fourier pseudospectral methodologies gives rise to a new methodology, named IMERSPEC, which was presented in [19]. The idea behind this method is that it combines the advantages of each technique and resolves some of their drawbacks. The base of this proposal and preliminary applications were given in [19], where the Taylor-Green vortex decay, lid-driven cavity flows and flows over square cylinders were solved using the IMERSPEC methodology. In [19], it was demonstrated that it is possible to solve non-periodic problems using Fourier spectral methodology, but the methodology is limited to simulate, only, flow over Cartesian geometries and presented low convergence rates-first order.
The listed works that have the purpose of evaluating some specific numerical and/or computational aspect in IMERSPEC methodology [20] are compared the analysis of convergence orders and computational times between the pseudo-spectral Fourier method and the finite volume method. In [21,22] is presented the development of the insertion of three different boundary conditions, using the immersed boundary method, particularly to thermal flow problems. In [23], the Lagrangian force is calculated by a function of the surface tension by the contact of two immiscible fluids, the immersed boundary is not rigid, and a different algorithm is developed. In the present work, the main goal is to show that the IMERSPEC methodology can be extended to more complex problems and achieve high-order numerical convergence. The results reach fourth order of spatial convergence and it is possible to establish that the convergence order can be modified by changing the distribution and interpolation functions from the immersed boundary method. Finally, applications to very complex flows over a backward-facing step are performed, and the two-dimensional and three-dimensional results agree with experimental results.

Mathematical modelling
In the present section, details of the merging of the immersed boundary [24], Fourier pseudo-spectral [8,13] and methodologies are presented. This coupling gives rise to the IMERSPEC methodology. The mathematical formulation is presented in the following subsections: mathematical modelling for the fluid flow, Navier-Stokes and continuity equations, both in physical space and in Fourier space, and interface equations, given by the immersed boundary and multi-direct forcing methods [25]. Finally, the coupling between the methodologies is presented.

Mathematical modelling for a fluid flow in physical space
In the present paper, mathematical models are presented for incompressible flows of Newtonian fluids without heat transfer and with constant physical properties. The formulation consists of the continuity equation, which arises from a balance of mass, and by the Navier-Stokes equations, which model the linear momentum balance (Newton's second law). These equations are presented using tensorial notation: In the present work, f i ( , t) is the force field of the fluid-solid interaction, the fluid-fluid interaction or another force field that models other physical effects internal to the fluid flow domain [3,4]. Equations 1 and 2, together with the boundary and initial conditions, can model any fluid flow in physical space, considering the simplifying assumptions used. In the present work, the Fourier pseudo-spectral method is used, and the mathematical modelling is presented in Fourier space.

Mathematic modelling for fluid flow in Fourier space
Defining the direct Fourier transforms, Eq. 3, and inverse Fourier transforms, Eq. 4, like [17,19]: where ̂ ( , t) are the Fourier coeficients of transformed function ( , t) , that can be the pressure or vector velocity, is the wave number vector, a parameter of the spatial Fourier transform and = √ −1. We make use of the definition of the forward Fourier transform (Eq. 3) in Eq. 1; it gives: where û( , t) are the coefficients of the direct Fourier transform of velocity vector. Equation 5 shows that velocity field for incompressible flows in Fourier space is orthogonal to the wave number vector. Application of the Fourier transform to Eq. 2 gives: where , and are the wave number vectors and k = ‖ ‖ is the norm of the wave number vector.
The term associated with the Fourier transform of the pressure field, k ip ( ,t) , is normal to velocity vector. The sum that composes the right-hand side of Eq. 6 is equal to the projection over plane divergent-free, of the sum of the nonlinear term (convolution integral) with the source term . Therefore, the transformed Navier-Stokes equations in Fourier space can be rewritten as: (5) k jûj = 0.
In Eq. 7 the projection tensor P im is defined as [26]: were im is the Kronecker delta tensor. We note that the P im tensor projects any vector over plane . It is noteworthy that, in Eq. 7, the pressure gradient no longer appears. This finding means that the velocity field, for incompressible flows, can be obtained without the pressure term. This property is one of the most interesting features of the Fourier spectral methodology. Therefore, with the proposed methodology, it is possible for complex problems to be solved without treating the pressure-velocity coupling and without solving linear systems.
Nevertheless, because the velocity field is solved, the pressure field can be determined, either as a direct result or as a post-processing result. The convolution integral, corresponding to the Fourier transform of a nonlinear term, which appears in Eq. 7, cannot be solved numerically. Thus, as an alternative, the pseudo-spectral method is described in the next section.

Fourier pseudo-spectral method
The Fourier pseudo-spectral method [7] involves the product of two functions in physical space and transforms this product rather than transforming the two functions independently and solves the convolution integral in Fourier space.
The advantage of this procedure is that we do not have to solve the convolution integral, although this reduces the accuracy of the Fourier spectral method. However, high accuracy can be maintained by performing the derivatives in spectral space and subsequently solving the product in physical space. The drawback of the pseudo-spectral method is that one must solve the direct and inverse Fourier transform for each time step, which increases the computational cost.
Another important detail regarding the presented formulation is the way that the nonlinear term is treated. The skew-symmetric form proposed by [7] is considered to be more stable and is used in the present paper. It consists of an arithmetic mean of the nonlinear term, written in conservative form and in non-conservative form.

DFT and FFT
Taking a discrete Fourier transform (DFT) is the proper numerical way to evaluate the forward and backward Fourier transform. However, it is restricted to periodic boundary condition problems, thereby limiting the use of numerical Fourier transforms for CFD problems. The Fourier spectral method has only been used for simulations of temporal jets, temporal mixing layers and homogeneous and isotropic turbulence [7,18]. The fast Fourier transform (FFT) algorithm proposed by [15] can be used to efficiently solve DFT equation. In terms of float point operations, the FFT gives O(N log 2 N) , whereas the DFT is on the order of O(N 2 ) float point operations where N is the number of collocation points. In the present paper, we use the FFTE subroutine implemented in [27].

Immersed boundary methodology-IBM
The IBM methodology uses two calculation domains, the Eulerian and Lagrangian domains. The fluid flow equations (Eq. 7) are solved in the Eulerian domain, and the force field is evaluated in the Lagrangian domain to represent the immersed interface of the flow [3]. These two domains communicate with each other through a source term added in the Navier-Stokes equations.
This source term is defined throughout the whole Eulerian domain with an understanding of its distribution over the Lagrangean domain, as given by: where F i ( k , t) is modelled according to that presented in [19]. The vectors and k are shown in Fig. 1.
The IBM used in the present paper was developed based on the work of [4], in which a balance of momentum is applied over an immersed interface (virtual physical model-VPM) by generating a reaction force that represents the boundary conditions of the immersed body. The IBM is also based on the work of [25], which reports an updated version of the VPM of [4]. This modification changes the manner in which the interface force is calculated through an iterative process, through multi-direct forcing (MDF), in order to improve the accuracy of the Lagrangian force calculus, F i ( k , t).
The domains used in the present work are presented in Fig. 1. The physical domain Ω PhD , termed the Eulerian domain, is bounded by the immersed boundary Γ PhD , which is considered as a Lagrangian domain. All physical phenomena that we are interested in modelling and simulating occur within Ω PhD . In Ω PhD , one or more bodies can be immersed, as represented by the frontiers Γ i , where i = 1, 2, … , N . Over the immersed frontier Γ PhD , any boundary conditions can be modelled even if they are not periodic, in the present article only Dirichlet boundary conditions are studied. Kinoshita et al. [21,22] show that thermal boundary conditions are assessed, [23] boundary conditions to two-phase flows are tested, and [18,28] turbulent wall models are presented.
The additional Eulerian domain Ω PeD is bounded by the frontier Γ PeD over which periodic boundary conditions are imposed. The domain Ω PeD complements Ω PhD and makes it possible to extend the solution for non-periodic boundary condition problems. Periodic boundary conditions are imposed directly over Γ PeD . On the other hand, non-periodic boundary conditions are imposed over the immersed boundary Γ PhD in an indirect way, through the Eulerian force field f i ( , t) given by Eq. 9.
In Fig. 1, the vector defines any point that belongs to the Eulerian domain, while the vector defines any point that belongs to the Lagrangian interface.
To evaluate the Lagrangian force F i k , t , we use the Euler scheme for temporal discretisation. From Eq. 7 we obtain: where t + Δt indicates the current time, t indicates the preceding time, u * i represents the i component of the predicted velocity, u t+Δt i and u t i are the solenoidal velocity (satisfying the mass balance), and rhs i is determined by summing the nonlinear, diffusive and pressure gradient terms. By decomposing Eq. 10 into two equations, the following EDP system is obtained: Equations 11 and 12 are valid for any fluid particle, including that positioned over the immersed boundaries. Next, Eq. 12 can be rewritten as: where F i is the Lagrangian force and U t+Δt i is the interface velocity, which is determined by the physical conditions related to the problem under study ("boundary conditions"). The velocity is zero for interfaces at rest, and it is nonzero for moving interfaces. The velocity for moving interfaces can be found by imposing a velocity for a rigid body, or it may not be known, e.g. in the case of fluid-structure interaction problems. In this case, the body velocity must be determined by the structural model.
The Lagrangian velocity U * i (predicted) is obtained by interpolating the predicted Eulerian velocity, u * i , obtained using Eq. 11.
On the other hand, the Eulerian force f i that appears in Eq. 12 is obtained from the distribution of F i , which is calculated using Eq. 13. Hence, the Eulerian velocity u t+Δt i can be updated using Eq. 14, which is given by Eq. 12: In order to improve the accuracy of the solution of the immersed interface velocity, the above procedure, i.e. the calculation of f i and u t+Δt i , can be performed in an iterative way. According to [25], the velocity field u t+Δt where it is one of the interactions of the multi-direct forcing method.
Next, u t+Δt,it+1 i is interpolated again, providing U it+1 i , and the Lagrangian force, F it+1 i , is recalculated and redistributed. Equation 14 is solved again, giving: This iterative process continues until a stopping criterion, , is satisfied. In the present paper, we adopted the force difference between two consecutive interactions: We note that, in each interaction, the velocity u t+Δt,it+1 i from Eq. 15 and the force f it+1 i are projected over plane , assuring mass and momentum balance. Furthermore, the force difference (Eq. 16) converges to zero, because the velocity at boundaries become more close to the boundary conditions. In Sect. 3.2.3, the results of this iterative process are presented and discussed.

Distribution and interpolation functions
It is noted that the force distribution and interpolation velocity processes, as described above, are given using the Dirac nucleus when the Lagrangian and Eulerian meshes are not coincident. In the present paper, fluid flows over Cartesian geometries are modelled. A wide range of problems in which the two meshes (Eulerian and Lagrangian) are coincident can be covered. In this case, the interpolation and distribution processes are replaced by the definitions given by Eqs. 17 and 18: However, when the problem has a non-Cartesian geometry, distribution and interpolation functions must be used. For this purpose, velocity interpolation and force distribution functions are defined. The following were proposed by [24]: where r x = x−x k h and r y = y−y k h , with h being the spacing between the collocation nodes of the Eulerian discretised domain, and Δs is the spacing between the nodes of dicretisation of the Lagrangian domain.
The W g function is similar in form to a Gaussian and satisfies the property of having a unit integral over [−∞, +∞] . The Eulerian force field, f i ( , t) , is null throughout the whole domain except in the neighborhood of the Lagrangian points where it virtually models an immersed interface, simulating the presence of a body. Thus, it is not necessary to make changes and adaptations to the Eulerian mesh in order to fit the mesh to the interface [4]. Once the Lagrangian force field, F i k , t , is calculated, it can be distributed and thus can transmit information about the geometry to the Eulerian mesh.
In the present work, in addition to the distribution function proposed by [24], i.e. Eq. 22, other functions proposed by different authors were also used. The main difference between the functions lies in the weight function, W(r). A cubic function, W c , was proposed by [29]; a six-point function, W 6p , was proposed by [30]; a cosine function, W cos , was proposed in the work of [24]; and a hat function, W h , was presented in [31].
Not all of these weight functions exhibit the mathematical properties that a distribution function must have. Nevertheless, in [32], it was shown that these functions give notably good results. The weight function, W(r), is fundamental to the accuracy and numerical convergence rates of the immersed boundary methodology.

Merging the pseudo-spectral and immersed boundary methods
The originality of the present work lies in one key idea: it uses the Fourier pseudo-spectral method to simulate non-periodic problems. This property is achieved by the IMERSPEC methodology proposed in [19], which allows for merging of the pseudo-spectral and immersed boundary methodologies. The algorithm used to solve the differential model is described below: 10. Determine whether is less than the specified tolerance, see Eq. 16. If the tolerance criterion is not satisfied, return to step 3. If the tolerance criterion is satisfied, go to the next time step (return to step 1).
We note that the Euler scheme, for temporal discretisation, was used only for demonstrating the algorithm. In the present work, the fourth-order, low-dispersion and low-storage Runge-Kutta method, which was developed by [33], is used. Furthermore, a variable time step was used based on CFL criterion [34].
where the CFL is a number between 0 and 1.

Manufactured solution
A synthesised or manufactured solution was used to determine a source term, thereby giving an analytical solution to the velocity and pressure fields. This source term, which has dimensions of [N/m 3 ] , is added to the Navier-Stokes equations, which has a numerical solution that can be compared with the synthesised solution [20]. The following equations, proposed by [35], are used in the present work: where u an and v an are the analytical solutions for the velocity components, p an is the analytical solution for the pressure field, a and b are constants, x and y are the components of (26) Fig. 2 Time distribution of the L 2 norm for both cases: the Taylor-Green problem with and without the immersed boundary the coordinate system, t is the time, and and are the fluid density and the kinematic viscosity, respectively. The main goal of this simulation was to verify the proposed methodology and its numerical implementation. In order to compare the numerical simulation with the analytical solution, the initial and boundary conditions were determined using Eqs. 27, 28 and 29. This procedure results in periodic boundary conditions, which are appropriate for the Fourier pseudo-spectral method. The control parameters, as suggested by [36], are the vortex diameters, L r = [m], and the maximum velocity, where Θ num and Θ an stand, respectively, for the numerical and analytical solutions of the generic variable field, Θ , i.e. (u, v, p).

Manufactured solution without and with a coincident immersed boundary
An initial simulation was carried out, setting the source term f m in Eq. 7 equal to zero. Figure 2 shows the time evolution of the L 2 norm for the u component of the velocity field. It shows the high accuracy of the Fourier pseudo-spectral methodology, as the maximum error reaches 10 −15 , though the simulation is performed with a coarse mesh of 16 × 16 collocation points. The second test for verifying the Fourier pseudo-spectral methodology was to place an immersed boundary in the Taylor-Green problem. This type of coupling was proposed and used in [36,37].
The first part of this study was done with the immersed boundary being coincident with part of the collocated Eulerian points, indicating the Γ PhD domain in Fig. 1 need in IMERSPEC methodology, a square of side length L∕L r = 1∕ is placed into Eulerian domain, Ω , and all of the Lagrangian points are coincident with the Eulerian collocation points. The procedure for verifying the coupling of the Eulerian and Lagrangian meshes involves setting the velocity of the immersed boundary U k , t + Δt equal to the analytical velocity calculated at the Lagrangian points.  The results of L 2 norm show that, in both cases, the maximum error is on the order of 10 −15 , as shown in Fig. 2, which characterises the high accuracy obtained with this methodology. It worth noting that the mesh used was 16 × 16 and that this calculation was also done for finer meshes with the same accuracy being obtained. This result characterises a methodology with spectral accuracy, as expected.

Manufactured solution with a non-coincident immersed boundary
This case was developed for a non-coincident immersed boundary to study the effects of the interpolation and distribution procedures given by Eqs. 19 and 20. A circular geometry of diameter L∕L r = 1∕ , where placed in centre of Ω domain and there is no coincidence between the Lagrangian points and the Eulerian points. The error must be attributed only to the interpolation and distribution procedures. Five distribution and interpolation functions were used, as described in [32], and the results show the behaviour of the methodology related to the type of the distribution and interpolation functions used. Figure 3 is obtained by doubling the number of Eulerian collocation points for each coordinate direction: 16 × 16 , 32 × 32 , 64 × 64 , 128 × 128 and 256 × 256 . The number of Lagrangian points was also changed: 28, 56, 112, 224 and 448 in order to maintain Δx ≃ Δs . For each simulation, the L 2 norm of the u velocity component is calculated after t∕t r = 1.0 . The time step Δt = 10 −4 [s] was used for all simulations. This time step does not influence the results obtained using different spatial refinements. In Table 1, the absolute maximum errors, of each function by the number of collocation points, are shown to the u velocity component.
The Gaussian, hat and cosine functions give the IMER-SPEC methodology a second-order convergence. Of these three functions, the hat function gives the most accurate results. The six-point function and the cubic function give the IMERSPEC methodology a fourth-order convergence and the cubic function gives the most accurate results. In Table 1 are shown the absolute maximum errors to five interpolation and distribution functions, where the high accuracy of the IMERSPEC methodology can be observed.
Importantly, all of the variables, including the pressure, show a numerical convergence of at least the fourth order. This is one of the main characteristics of the IMERSPEC methodology: it is an immersed-boundary-based methodology with spatial convergence of the fourth order. Another important characteristic is that the pressure field presents the same order of convergence as the velocity. Normally,  conventional methods present an order of convergence for the pressure that is smaller than that of the other variables. The result for the velocity v presents the same error and the same order of convergence as the velocity u. It is interesting that the Lagrangian velocity U also presents a convergence of the fourth order. This result is important and shows that the flow is not affected by the presence of the Lagrangian immersed boundary.
Comparing the presented results with the reports published up to this point, the proposed methodology gives one of the highest orders of convergence when the immersed boundary method is used. The works of [36,37] presented simulations of the same problem with the same parameters. They obtained second-order convergence and a low level of accuracy, as compared with the present work's results obtained using the interpolation/distribution methods given the six-point function and the cubic function.
However, it is important to note that manufactured solutions used (Eqs. 27, 28 and 29) are smooth and the Gibbs phenomenon does not appear, even in simulations that immersed boundary is used, and there is no coincidence between Lagrangian and Eulerian points. The present results characterise that IMERSPEC methodology has been verified. When there is coincident points, it reaches round-off errors, and when there is not coincident points, the fourth order of spatial convergence is obtained. The limitations of this order of convergence are the interpolation/distribution functions and the temporal integration method. The next test case is not smooth, and the Gibbs phenomenon will be again addressed.

Two-dimensional backward-facing step flows
The backward-facing step is one of the most commonly used test cases for verifying the potential of a new methodology. It is geometrically simple and presents a high control of the boundary layer detachment point. On the other hand, it generates a notably complex flow with detachment and reattachment points, shear flow, recirculation zones and boundary layers. Many related experimental reports have been published, for instance, [38][39][40], as well as numerical results [41][42][43][44]. Figures 4 and 5 show, respectively,the geometrical information and the expected physical characteristics of the flow over this geometry. The geometry is characterised by an entrance channel of length L in and height (W − h) . The step has a height h. The channel after the step has length L and height W. The ratio W/h is known as the aspect ratio. Figure 5 illustrates the boundary layers at the inlet channel. The detachment of the boundary layer at the bottom wall occurs at (x = 0, y = h) independently of the flow upwind of this point. The boundary layer reattaches at (x = X r , y = 0) . The boundary layer at the top wall detaches at (x = X s , y = W) and reattaches at (x = X s + X rs , y = W).
The detachment of the bottom wall boundary layer promotes a recirculation zone and creates a mean inflectional velocity field, which gives rise to Kelvin-Helmholtz instabilities. The adverse pressure gradient downstream of the step causes the top wall boundary layer to detach and creates Kelvin-Helmholtz instabilities below the top wall. Therefore, the main characteristic of this problem is that the geometry is simple, but the flow inside it is complex, and within it one can find several kinds of physical instabilities, such as boundary layer instabilities, Kelvin-Helmholtz instabilities, collisions of instabilities with the walls, interactions between the instabilities, boundary layer detachments, boundary layer reattachments and boundary layer and Kelvin-Helmholtz interactions. The interactions of Kelvin-Helmholtz instabilities with the walls create counter-rotating pairs that can cross the entire channel, going from one wall to the opposite one. Therefore, this geometry results in a complex and interesting benchmark for validating a new methodology. Figures 4 and 5 show that this problem is not periodic in any direction; the main objective of the present work is to make it periodic using the IMERSPEC methodology.
This goal is the natural step following the work of [19]. The procedure consists of enclosing the non-periodic physical domain within a periodic domain, which is done by using the immersed boundary methodology to force the real nonperiodic boundary conditions inside the periodic domain. Figure 6 shows the periodic domain Ω PeD , which is delimited by the periodic boundary Γ PeD , which is represented by the dashed line. The physical domain Ω PhD , which is delimited by the physical boundary Γ PhD , is represented by the solid line.
Over the entire frontier Γ PeD , periodic boundary conditions are directly imposed. Over the frontier Γ PhD , one can impose any kind of "boundary conditions" using the immersed boundary method.   "boundary condition" is imposed, representing the presence of solid walls. From the point (0, h) to the point (0, W), a known velocity profile is forced, including the residual turbulence, which is measured experimentally. Importantly, this alternative way of imposing a "boundary condition" uses the source term from the Navier-Stokes equations. This is a Eulerian force term, ( , t) , which is obtained by the distribution of the Lagrangian force field, k , t . As previously described, k , t is calculated by the multi-direct forcing method.
For the present work, special care was taken to obtain convergence and physical consistency. Due to its periodicity, the physical instabilities that pass through the outlet of the domain are re-injected at its entrance, as illustrated in Fig. 7.
These instabilities must be dampened using a buffer zone BZ, as illustrated in Fig. 6. The dampening process ensures that the target velocity is attained at the entrance of the forcing zone FZ. The forcing zone FZ is used to align the stream lines before entering the physical domain Ω PhD . This process is enforced in every case, assuring that the target velocity is imposed over the entire domain FZ. In the present work, the mean velocity profile to be forced is as follows: This equation, which models the mean velocity at the physical domain entrance, is the same as that proposed in [38], which corresponds to a developed flow inside a rectangular channel. The mean velocity is Ū = 1.0 [m/s]. The dampening process in the buffer zone BZ is carried out using the dampening function given by [45]: where Q t i is the i component of the target velocity. In the present case, one takes Q t x = U in (y) in the streamwise direction and Q t y = 0 in the vertical direction. Q i is the solution obtained from the Navier-Stokes equations. The function , which models the vortex strain, is given by: where x sB and x fB are, respectively, the beginning and the end of the buffer zone.
All of the simulations presented in the present work were carried out using the Runge-Kutta method, which is of fourth order, using the optimised coefficients proposed by [33]. A value of CFL = 0.5 (Eq. 26) was also used. The simulation domain, presented in Fig. 6, has the following dimensions: h = 0.5 [m]; L x ∕h = 73.14 ; L y ∕h = 2.29 ; W∕h = 2.0 ; L BZ ∕h = 3.73 ; and L FZ ∕h = 0.53 . The kinematic viscosity is given by = hŪ∕ℜ h , where ℜ h is the Reynolds number and Ū is the mean velocity at the entrance of the physical domain.

Two-dimensional backward-facing step at ℜ h = 400
In the present section, the simulated flow refers to ℜ h = 400 , and the previously defined domain is divided into N x = 2048 and N y = 64 collocation points. In defining control parameter for the multi-direct forcing method, a maximum cycle number of NL = 200 or a maximal residual (Eq. 16) of = 10 −5 were used. Figure 8a shows the streamwise velocity component, and Fig. 8b gives a detailed view of streamlines.
The classical laminar flow with the bottom and the upper recirculation zones is visualised in Fig. 8a, b. These recirculations appear because of the detachment of the boundary layers. The laminar behaviour for this Reynolds number is as expected. To demonstrate the importance of the buffer and forcing zones, Fig. 9 shows the manner in which the fluid particles are driven to the entrance of the physical domain Ω PhD . Figure 9a shows the streamlines for the case with the forcing zone where we can see that they are parallel to the walls, as must happen in a wind tunnel, as illustrated in [19]. Figure 9b shows the streamwise velocity component values at y∕h = 1.49 obtained with and without the forcing zone. This figure also shows a reference solution or the target velocity, which is given by Eq. 31. We see that the streamlines become parallel for the fluid entry into the forcing zone with the target velocity field used as a reference to calculate the force profile given by Eq. 31. Without the forcing zone, the velocity at the entrance plane does not reach the target velocity value and is not parallel to the walls. The oscillation that appears after x∕h = 0.0 is already inside the domain   and is the result of the first physical instabilities, such as the Kelvin-Helmholtz instabilities. Figure 10 presents comparisons between the results of the present work and the experimental results of [38]. The streamwise velocity is shown at three points: (a) x∕h = 0.0 ; (b) x∕h = 14.0 ; (c) x∕h = 30.0.
At point (a), very good agreement is obtained, showing that the internal boundary is well modelled by the forcing zone. At points (b) and (c), again, very good agreement is attained with the experimental results of [38] and with the numerical results of [41]. The L 2 norm was also calculated at the entrance of the domain that measures the difference between the numerically recovered velocity and the target velocity, which is given by Eq. 31, at the entrance of the backward-facing step. The maximum value obtained for this norm was L 2 = 10 −5 [m/s] for both velocity components, as shown in Fig. 11. This result shows that the method used to force the internal "boundary condition" is notably accurate and stable during the time evolution of the numerical experiment. It is possible to accurately reproduce the experimental measurements, regarding the numerical representation of the boundary conditions.

The influence of the mesh refinement
The next simulations and results show the behaviour of the solution as the mesh is refined. The same domain and the  same parameters for the multi-direct forcing method used in the previous section are used in the present section. Three different simulations were carried out using three different levels of mesh refinement. Importantly, the L 2 norm measures how close the forced boundary condition is to the proposed boundary condition. For instance, the L 2 norm that is calculated at the entrance measures how close the forced boundary condition is to the experimental velocity profile.
On the other hand, the L 2 norm calculated at the horizontal walls measures how close we are to the no-slip boundary condition. Figure 11 show that the forcing process is sensitive to the mesh refinement and that the accuracy of the method increases as the mesh is refined.
The skin friction coefficient, C f , is a important parameter measure in backward-facing step [46][47][48], and is expressed by Eq. 34: where w is the wall shear stress, w = du dy y=w In Fig. 12 shows the distribution of skin friction coefficient downstream of a backward-facing step on bottom (a) and top (b) walls, for the three different meshes. The results of Fig. 12 are used to obtain the position of reattachment point for the bottom wall, and the separation and reattachment points for the upper wall, respectively, are given in Table 2. Table 2 presents the reattachment point for the bottom wall and the separation and reattachment points for the upper wall. These parameters are compared with the experimental results of [38] and the numerical results of [41]. The convergence of results is highlighted, i.e. when the refinement increase the results close to experimental and numerical results. The results of the present work improve as the mesh is refined. The approximate CPU time and the RAM memory used in each simulations are provide in Table 3; these values are obtained on a computer Intel ® Xeon® CPU E3-1270 v3@, 3.5 [GHz] and 16 [Gb] RAM memory. The operational system is the UBUNTU 14.04. The result of CPU time shows that the relationship of the total collocation points is approximately K ⋅ N ⋅ log 2 N , where K is a constant and N is the total of collocation points.

The multi-direct forcing method
In this section, the influence of the number of loops (NL) on the multi-direct forcing method is presented. A discretisation process is used with N x = 2048 and N y = 64 collocation points. The parameter NL was used as the simulation stopping criterion. Figure 13 shows the "internal boundary condition" at the domain entrance, x∕h = 0.0 (see Fig. 6), and at the first time step, t = Δt . The velocity profiles are shown for different interaction numbers using the same simulation. As the NL parameter is increased, the forced velocity profile becomes independent of this parameter and tends to the experimental velocity profile measured at this section of the domain, as given by Eq. 31.
The difference between the forced and the experimental inlet velocity is evaluated using the L 2 norm calculated using all of the collocation points placed at the entrance plane. Figure 14a shows the time distribution of L 2 for the u component of the velocity at the inlet of the domain. Similarly, Fig. 14b shows the L 2 norm at the bottom wall. The difference between the forced velocity and the target velocity is slight, and, as the cycle number of the direct forcing is increased, the L 2 norm reaches 10 −5 [m/s] and 10 −7 [m/s]. Thus, the "internal boundary conditions" are well represented when compared with the experimental results at the entrance and with the no-slip condition at the solid walls. Table 4 shows the influence of the cycle number of the direct forcing method for the following parameters: the reattachment point at the bottom wall, x r ∕h ; the separation point at the upper wall, x s ∕h ; and the separation point at the upper wall, x rs ∕h . It also presents a comparison of the experimental and numerical results of [38,41]. We see that the cycling process is notably important. With only one cycle, the error is considerable. Increasing the cycle number decreases the error, which reaches a notably low level. The results show that, for NL greater than 10, the results do not change significantly.

Different Reynolds numbers
In this section, the simulation results for several Reynolds number values are presented. A mesh of 2048 × 64 collocation points was used. The Reynolds numbers studied range from 200 to 6000. Figure 15 presents the vorticity fields for the backward-facing step. The buffer and forcing zones are not presented. For Reynolds numbers ranging from 200 to 600, the flow does not present any instability, as expected. For Reynolds numbers ranging from 800 to 2500, the flow becomes instable, giving rise to Kelvin-Helmholtz instabilities near the step that are transported towards the bottom wall reattachment point. Instabilities are also formed beyond the detachment point on the upper wall. These instabilities interact with each other and are transported towards the outlet of the channel. This finding is a classical result that is obtained with the proposed methodology in the present work. For Reynolds numbers ranging from 3000 to 6000, the instabilities present chaotic behaviour, as expected for higher Reynolds numbers. Table 5 shows the reattachment point at the bottom wall ( x r ∕h ), for several Reynolds number values and a comparison with the experimental results of [38]. There is good agreement for Reynolds numbers smaller than 400.
For higher values of this parameter, the error increases as expected for two-dimensional simulations, as was also reported in [39]. Figure 16 shows the details of the flow, which gives rise to considerable information. Just downstream of the step, the shear layer is captured, which gives rise to the Kelvin-Helmholtz instabilities, as shown in Fig. 16 (a).
Afterwards, there are strong interactions between the "turbulent" structures and the walls, as shown in Fig. 16b. Notably interesting physical effects are observed in Fig. 16c. Figure 16d shows the details of the vortex injection inside the domain caused by very strong shear effects. These types of physical effects and details can only be obtained when highly accurate numerical methods are used.
Furthermore, the instabilities that leave the domain, Fig. 6, are damped by buffer zone and the function , (Eq. 33) is fundamental to smooth the outlet flow and avoid the Gibbs phenomenon.
We would like to stress that no turbulence model was used; thus, only methods with high accuracy or mesh refinement can capture the flow details presented in Fig. 16.

Three-dimensional backward-facing step flows
Furthermore, the IMERSPEC methodology has been extended to study the three-dimensional fluid flows and was developed a computational parallel code, which is verified in work of [18]. This three-dimensional computational code has the same features of the two-dimensional computational code and presents three turbulence models based in large eddy simulation (LES) and it is parallised by using distributed memory multi-computer with MPI. For more details, see the work of [18]. In the present work, the three-dimensional parallel code is used to analise the backward-facing step flows.
The simulation domain, presented in Fig. 6, has the following dimensions: h = 0.5 [m], L x ∕h = 54.86 , L y ∕h = 2.29 , W∕h = 2.0 ; L BZ ∕h = 3.60 e L FZ ∕h = 0.70 , with NL = 200 number of loops on the multi-direct forcing method and the value of CFL = 0.5 (Eq. 26) was also used. The third dimension has length L z ∕h = 9.16 . The inlet velocity profile, given by Eq. 31, is imposed throughout z direction in forcing zone ( L FZ ∕h ) and periodic boundary conditions are directly imposed in z direction. The domain is discretised using N x = 768 , N y = 32 and N z = 32 collocation points and to the parallel computation the domain is divided equally by four processor in x direction, i.e. each processor gets N x = 192 , N y = 32 and N z = 32 collocation points.
In this work simulated the three-dimensional backwardfacing step flow at two Reynolds numbers, ℜ h = 400 and ℜ h = 1000 . In the simulation at ℜ h = 400 , a comparison of horizontal velocity profiles of the two-dimensional results  of the present work (IMERSPEC 2D) and the experimental work of [38] is presented in Fig. 17. In Fig. 17 at points (a) and (b), very good agreement is attained with both the experimental results of [38] and the two-dimensional numerical results of the present work. The results of the velocity in 2D and 3D simulations are close, as expected, due at ℜ h = 400 is laminar flow and there is not velocity component in z direction, i.e. w ≈ 0.
In Fig. 18 is presented the plot of vorticity field, z , in different viewing positions to the simulation of three-dimensional backward-facing step flow at ℜ h = 1000. Figure 18 shows that a recirculation appears because exist the detachment of the boundary layer. From the reattachment point in the bottom wall arises the transition to turbulence flow. Notably, interesting physical effects can be observed, which gives rise to the Kelvin-Helmholtz instabilities, as shown in Fig. 18a, c. The Kelvin-Helmholtz instabilities cause the formation of eddies, which interact with the botton and upper walls. In Fig. 18b-upper view, the vortex begins to oscillate in z direction, one of the characteristics of the transition to turbulence. Furthermore, Table 6 shows the comparison of three works of the mean position of the reattachment point for the bottom wall ( x r ). This parameter is compared with the experimental results of [38] and the two-dimensional numerical results of the present work (IMERSPEC 2D).
The simulation at ℜ h = 400 to IMERSPEC 3D presents the position of the reattachment point for the bottom wall, x r , close to the results of the simulations of the IMERSPEC 2D and experimental data of [38] as expected, because the flow is laminar and stable. The three-dimensional numerical results at ℜ h = 1000 are close to, only, the experimental data [38]. In two-dimensional simulations, even using a refined mesh ( 2048 × 64 collocation nodes) the results do not come close the experimental data. These results show that the backward-facing step flow at ℜ h = 1000 has strong three dimensional characteristics. Then, the three-dimensional simulations is required to analise this flow.
Lastly, Fig. 19 shows the time distribution of L 2 norm for the horizontal velocity component at the inlet of the domain (a) and at the bottom wall (b) to different Reynolds numbers.
The errors present in Fig. 19 are larger in the inlet plane, order of 10 −4 , due to the influence of the vortex that appears because the periodic bondary coditions, to see Fig. 7. In the botton wall, the errors are smaller order ( 10 −5 ) and the errors of the upper wall and of the others velocity components (v and w) were not shown, but are the same order of the errors of the botton wall.

Conclusions
In the present paper, the IMERSPEC methodology was presented, which is based on a merging of the immersed boundary and Fourier pseudo-spectral method, showing to be possible to use Fourier pseudo-spectral method to simulate flow over complex geometries retaining part of its accuracy. The immersed boundary method is based on the works of [4,25]. This methodology may be applied to a large number of problems that can be described by a Cartesian immersed boundary. For non-Cartesian problems, fourth-order convergence is attained. The proposed methodology was rigorously verified through comparison with a synthesised analytical solution. Furthermore, a weighted spectral accuracy was attained for a large spectrum of Cartesian problems. Thus, after the numerical verification, a validation work was carried out.
The complex flows occurring over a backward-facing step were chosen. Two-dimensional simulations were carried out, and the results were compared with the experimental and numerical results of other authors. Very good agreement was obtained, and we consider the proposed methodology to be promising for solving complex flows.
The more important characteristics of the proposed methodology include that the mass balance is on the order of the round-off errors; that there is no velocity and pressure coupling for incompressible flows, and as a consequence, that no linear system needs to be solved for the pressure; that the pressure field can be recovered as a post-processing procedure, giving the same accuracy and the same order of convergence as the velocity field; and that the methodology is easy to use and implement.