Caroline Collange

David Defour
email: david.defour@univ-perp.fr

Arnaud Tisserand
email: arnaud.tisserand@irisa.fr

Power Consumption of GPUs from a Software Perspective

GPUs are now considered as serious challengers for highperformance computing solutions. They have power consumptions up to 300 W. This may lead to power supply and thermal dissipation problems in computing centers. In this article we investigate, using measurements, how and where modern GPUs are using energy during various computations in a CUDA environment.

Introduction

As GPUs gained in flexibility through high-level languages such as CUDA, GPUs gained interest for the acceleration of tasks usually performed by a CPU thanks to the high computational power of GPUs. Therefore we are witnessing a tremendous growth in the usage of GPUs for high-performance solutions in computing centers. However, as long as the main goal of GPU was to serve in desktops, their power consumption was secondary. Even though that the ratio of Single Precision GFLOP per W is always in favor of GPUs compared to traditional CPUs (4 SP GFLOP/W for a GTX280 and 0.8 SP GFLOP/W for a core i7 960), these processors are known to have high power consumption. Therefore new challenges need to be solved in order to spread their usage in computing centers where 1 dollar spent in power supply corresponds to 1 dollar spent in cooling system. This work investigates how and where the power consumption is located within a GPU board by analyzing the relations between the measured power consumption, the required time and the type of units that are stressed to perform a defined operation. There are several solutions to measure or estimate the power consumption of a processor. There exist power estimations based on cyclelevel simulation like in Wattch [START_REF] Brooks | Wattch: a framework for architectural-level power analysis and optimizations[END_REF] or SimplePower [START_REF] Ye | The design and use of simplepower: a cycle-accurate energy estimation tool[END_REF] that rely on a low-level description of the architecture which is unavailable for current GPUs. Functional approaches and tools such as SoftExplorer [START_REF] Senn | SoftExplorer: Estimating and optimizing the power and energy consumption of a C program for DSP applications[END_REF] to estimate power consumption at assembly or C levels were proposed for VLIW processors or DSPs. Another solution is to physically measure the power consumption. The measurements can also be used to calibrate some models used for power estimation. Previous work have shown that in the context of multimedia applications for modern out of order processors, CMP is more energy efficient than SMT [START_REF] Sasanka | The energy efficiency of cmp vs. smt for multimedia workloads[END_REF]. [START_REF] Weste | CMOS VLSI Design: A Circuits and Systems Perspective[END_REF] and [START_REF] Piguet | Low-Power Electronics Design[END_REF] provide good references for power consumption and reduction aspects.

In this paper we consider Nvidia GPUs used for GPGPU (General Purpose computing using GPU) in a CUDA environment. During the analysis, functional blocks are identified, and their consumption is characterized using physical measurements. The considered blocks correspond to units that are usually stressed while executing common kernels on the GPU: register file, memory hierarchy and functional units. In addition to the power estimation, our analysis gives us some information on the organization of the memory hierarchy, the behavior of functional units and some undocumented features.

A modern GPU board is described in Section 2. The measurement process and the parameters extraction method from CUDA are described in Section 3. Results and their analysis for several GPUs are presented in Section 4.

Description of a GPU card

High-end GPUs considered in this work are sold with other component on a dedicated graphics card. These components consist in the GPU, the graphics BIOS, graphics memory, digital to analog converters, dedicated accelerators, an interface with the motherboard, a cooling device and power adapters.

While CPU manufacturers have moved toward higher power efficiency, the power required by GPUs has continued to rise. Although the power provided by power supplies have followed the necessary trend, a PCI-e add-in card can draw a maximum of 75 W through the standard connector, as specified in PCI Express CEM 1.1. GPU manufacturers provide extra power to the GPU by now using 8-pin (6-pin) wire-to-board connector that provide up to 150 W (75 W) additional power. This leads up to 300 W for one GPU card.

Description of Nvidia GPUs

In this paper we consider GPUs that are compatible with the DirectX 10 standard which provides a unified architecture where vertex and pixel shaders share the same instruction set and processing units. In order to efficiently exploit data parallelism, GPUs include a large number of replicated copies of these units operating in a single-instruction multiple-data (SIMD) fashion. GPUs handle high-latency instructions such as memory accesses by overlapping them using thousands of threads.

The unified architecture proposed by the DirectX 10 standard has been implemented in hardware since the release of the NVIDIA GeForce 8. An example of this architecture is depicted in Figure 1. The GPU has its own memory which bandwidth is usually an order of magnitude higher than system memory bandwidth. The graphics processor is seen as a set of multiprocessors. Each multiprocessor consists in numerous processing elements (PEs) operating in SIMD. At each clock cycle, all the PEs in a multiprocessor execute the same instruction sequence, but operate on different data. These multiprocessors incorporate various types of memory such as a register file, a scratchpad memory shared by all the PEs in a SIMD block, and read-only constant cache and texture cache.

In addition, PEs can also read or write global memory available on the graphic card.

The GeForce 8 has a shared cached constant memory accessible in broadcast mode only, a global memory accessible with coalescing (with additional alignment constraints), and a shared memory that allows coalesced, broadcast, and other patterns. Memory accesses that do not match these patterns are replaced by as many serial accesses as necessary, resulting in decreased performance. For example, the Tesla C870 embeds 1.5 GB of global memory with a peak bandwidth of 77 GB/s, which is more than 10 times the bandwidth available between a CPU and system memory. Multiprocessors integrate various computational units in order to implement the functionalities offered by the shaders: general computational units with embedded multiply-accumulators, texturing and filtering units, a dedicated unit to evaluate mathematical functions (e.g., sine, cosine, reciprocal, and reciprocal square root), and attribute interpolation units. Computational units can handle integer and floating-point arithmetic and there is no overhead associated with mixing both operations. Each multiprocessor of the GeForce 8 is able to execute a warp of 32 floating-point additions, multiplications, multiply-and-adds or integer additions, bitwise operations, comparisons, minimum or maximum of 2 numbers in 4 clock cycles. As there is no 32-bit integer multiplication in hardware, evaluating such an operation requires 16 clock cycles for a warp. For instance, one GPU of the GeForce 9800GX2, depicted in Figure 1, has 16 multiprocessors and each SIMD block is composed of 8 PEs and 2 functional evaluation units.

Measurement and Tests Description

We measure the power consumption of several Nvidia graphics card supporting the CUDA API described in Table 1. This includes a Tesla C870 card used for scientific computing, a dual GPU high-end graphics card 9800 GX2 and a T10P compute processor engineering sample. Table 1. Main characteristics of the tested graphic cards. Frequencies of the GPU, computational units and memory are given, as well as the manufacturing process, the number of transistor and the temperature provided by the integrated sensor.

The considered GPUs are tested in a desktop environment with a MSI X48 Platinum motherboard, an Intel Core 2 Duo E8400 processor with 6 MB L2 cache running at 3 GHz, four 1GB Corsair TWIN3X2048-1333C9DHX DDR3-1333 9-9-9-24 memory modules, and a 750 W Corsair HX750W power supply. This configuration is described in Figure 2. Our motherboard based on an X48 Intel chipset offers two PCI-e x16 2.0 and two PCI-e x4 1.0 slots through a bridge. The tested GPUs are included in a separated Tesla D870 box connected to the system through a dedicated bridge (NForce 100). Our software configuration uses Linux Ubuntu 8.04, with CUDA 2.0 and proprietary graphic drivers Nvidia 177.13.

To measure the power consumption we use a Tektronix TDS 3032 oscilloscope with the built-in analogue 20 MHz low-pass filter. The first input measures the current with a clamp sensor CA60, while the second measures voltage. The product of both measurements gives us the power consumption of the GPU through the external power. We took into consideration the power provided to the PCIe bus as well as the external power of the GPU. We noticed during our tests that the main source of power is provided by the external link, as during computation only 10 to 15 extra W are coming from the PCI bus.

We are aware that the measurements collected using this methodology are subject to caution. First, we not only measure the power consumption of the GPU but we also measure the power consumption of the whole board that includes the GPU, the graphic memory, DC/DC voltage converters, and others ICs (bridge, video chipset, . . . ) However we assume that we will notice a difference in the power consumption only if these parts are stressed. Secondly, we measure the consumption ahead of the board voltage regulators which include Fig. 2. Overview of our test configuration. small capacitors (decoupling). This means that the power consumption may be missestimated. We design short tests such that some variations in ambient temperature or the cooling system will not affect the results. Sampling rate is 50 kHz which allows us to measure the power consumption at task level but not at instruction level.

Results

Global Power Consumption

We measure the computation time and the average power required by common GPGPU algorithms (matrix transposition, matrix multiplication and cuBLAS) on 1024 × 1024 random matrices. The implementations of GPGPU algorithms used for these tests are the ones included in the Nvidia CUDA SDK 1.2. Results are reported in Table 2.

We observed during our tests that before the Linux driver was loaded, the GT200 was using around 20 W, which corresponds to the power saving mode, as officially claimed by Nvidia. Once the driver is loaded, the power consumption of the GT 200 rises to 51 W while idle. Surprisingly, we noticed that even though the GT200 includes low-power GDDR3 compared to the Tesla C870 and a smaller fabrication process, the GT200 requires more energy than the Tesla C870 to perform a naive or an optimized matrix transposition. As the energy depends on memory access patterns on data which are located in memory bank modulo the bus width, the example choosen (matrix of size 1024) is in favor of the C870 with a bus width of 384 bits compared to the 512 bits of the GT200 that raise bank conflicts when accessing data. Nevertheless, the GT200 requires twice less energy than the G80 on the matrix multiplication or cuBLAS example. This suggests that memory access patterns play a major role from performance as well as power consumption perspectives.

Multiprocessor

Considered GPUs integrate several multiprocessors with their own front-end: instruction fetch, decode, issue logic and execution units, etc. We run tests where the number of active multiprocessors varies from 1 to the maximum (16 for G80 and G92 and 30 for GT200). Results are reported in Figure 3. We observe that power consumption rise linearly with the number of multiprocessors up to respectively the maximum, one half and one third of the available multiprocessor for the G80, G92, GT200. This implies that the scheduling strategy for blocks to execute is different between the G80 and the G92 / GT200. We can deduce that the scheduling strategy for G92 and GT200 for blocks execution is to execute one block per group of texture processors which includes 2 multiprocessors for the G92 and 3 multiprocessors for the GT200. This means that the scheduling strategy of the G92 and GT200 is optimized for bandwidth usage, whereas the one used in the G80 is optimized for energy saving when the executed kernel is computationally bounded.

Execution Units

We measure the power consumption of various units within the GPU. We compare the power consumption of instructions predicated to true (that are really executed) and instructions that are predicated to false and we notice a significant variation. This means that instructions predicated to false are fetched, decoded but are not executed. This is unlike scalar architectures such as ARM or IA-64 behave, were predication is used to avoid pipeline stalls. On GPUs, predication is only used for SIMD control flow, while pipeline stalls are avoided using multithreading. We measure the power used with various combination of units (MAD units, MAD and MUL units, MAD and complex function evaluation units and double-precision unit). Results are reported in Table 4.

Memory Hierarchy

We measure how the memory hierarchy impacts the power consumption of GPUs. To measure the consumption, we developed a pointer chasing benchmark to access each level of the memory hierarchy. This benchmark includes a loop that traverses an integer array with various strides. 3 consecutive reads of one block of 128 bytes that will exploit coalescing are measured. Time is measured using the internal cycle counter of one multiprocessor. Only one warp on one multiprocessor is used. The size of the array and the stride vary. In order to minimize measurement noise, the minimum values among 10 executions are reported. Results for memory latency are given in figure 4. Table 3. Energy per memory request depending on where data is located. Energy results (nJ) per memory request of 128 bytes correspond to the measured power consumption divided by the measured bandwidth.

256 1K 4K 16K 64K 256K 1MB 1K 16K 256K 4MB
The observed latency for each GPU varies between 300 and 800 ns. In Figure 4.a (C870), we observe 3 steps at respectively 345 ns, 372 ns and 530 ns depending on the locality of the access. G92 behaves the same way with respective latencies of 325, 350 and 500 ns. We believe that these variations come from TLB hits and misses. Therefore, these results suggest a size of 4 KB per page with a 32-entry fully associative TLB. Results suggest that the second TLB level is accessed in 16 cycles on the G80 and G92 and 18 cycles on the GT200. However, our tests do not allow us to determine if it consists in an actual TLB or of a page table entry (PTE) cache. It is possible that it is shared with the second level of constant cache included in each texture processor. This cache may be used for instructions, constants, and PTEs.

DRAM is divided into pages that need to be precharged before each access and subsequently unloaded. According to our tests, DRAM page misses do not impact measured latency and power consumption. This may come from the fact that the memory controller optimizes page activation and unloading according to the address and type of memory access waiting in the buffer and that the access latency is made constant based on the worst case to simplify the design of the memory controller.

The latencies of texture memories are given in Figure 5. We observe the same steps as in Figure 4 plus the variation due to the two texture caches. By comparing the latency of a cache miss in the texture cache with a memory read in the same conditions, we estimate at 30 to 50 cycles the overhead due to texture filtering and the extra cache access depending on the GPU generation. We also measured that the latency to access the L2 cache and to go through the memory crossbar should be between 35 and 50 cycles. We measured the energy required for one memory request of 128 bytes executed by a warp on the G80, G92 and GT200 architectures. Results are reported in Table 3. The G80, G92 and GT200 are respectively made of 16, 16 and 30 multiprocessors. Kernels used for the tests are launched with 512 threads/block with as many blocks as multiprocessors. The instruction count reflects the number of operations that can be executed in a single warp. The measured power corresponds to the average power measured during the whole computation. CPI gives the number of clock cycles per instruction which corresponds to the number of clock cycles measured with the internal cycle counter divided by the number of threads per block. The energy per warp corresponds to extra power required to perform a specific operation for a warp multiplied by the computational time. We observe that the L1 texture cache of the G80 requires the same amount of energy per byte as the L2 texture cache. This is certainly due to the small bandwidth performance of the L1 texture cache (76 GB/s for L1 and L2 texture cache). However we cannot explain the small performance of the L2 texture cache of the GT200. The comparison of Table 3 and 4 leads to the conclusion that executing MAD is 7 to 15 times more energy efficient than accessing memory on a G80.

Conclusion

This article presents some measurements and an analysis on how computations and memory accesses impact the power consumption of some Nvidia GPUs (G80, G92, GT200) in the CUDA framework. Memory accesses can significantly degrade the performance and power consumption such that the G80 may be a better alternative than the latest GT200 in some specific cases. Our tests show how blocks of threads are dispatched among the multiprocessors when the number of multiprocessors is larger than the number of blocks and how it could negatively impact the power consumption for computationally bounded kernels. We also compare the energy required for various operations for a given warp with the energy required for a memory access and show that computations are up to 7 to 15 times more power efficient than memory accesses.

Our tests show that environmental conditions such as temperature and measurement equipment can impact the results by 10 to 15%. Therefore we are currently working on minimizing the impact of the environment by executing the same tests with others tools in a different environment. We also plan to measure the clock gating impact and perform comparisons between Nvidia and ATI GPUs for a given workload.

Fig. 1 .

 1 Fig. 1. Diagram of G92 architecture.

Fig. 3 .

 3 Fig. 3. Power variations depending on the number of active multiprocessors.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Latency for one read in global memory

Table 2 .

 2 Average power consumption P Avg , computation time T and corresponding energy E measured for the execution of common GPGPU algorithms (Naive and optimized matrix transposition, matrix multiplication, cuBLAS).

	GPU Idle Trans. naive Trans. optimized	MatMul	cuBLAS sgemm
		PAvg T	E PAvg T	E PAvg T	E PAvg T	E
		(W) (W) (ms) (J) (W) (ms) (J) (W) (ms) (J) (W) (ms) (J)
	G80	68 103 2.4 0.247 127 0.30 0.038 132 25 3.3 135 11.8 1.59
	G92	71 94 3.66 0.344 105 0.45 0.047 117 23.3 2.73 122 11.4 1.39
	GT200 51 73 4.11 0.300 83 0.50 0.042 114 13 1.48 113 7.44 0.84

Table 4 .

 4 Number of instructions, measured power P, number of cycles per instruction and the energy required per warp for ALU operation on G80, G92 and GT200.

			G80	G92	GT200
	Operation	#	P CPI	E	P CPI	E	P CPI	E
		inst. (W)	(nJ/warp) (W)	(nJ/warp) (W)	(nJ/warp)
	MAD	32 107 4.75	8.57	100 4.29	5.06	91 4.3	5.31
	Pred	32 90 2.38	2.43	93 2.39	2.14	75 2.36	1.75
	MAD+MUL 64 117 3.19	7.24	111 2.83	4.61	102 2.82	4.44
	MAD+RCP 40 115 3.96	8.63	110 3.55	5.63	98 3.54	5.14
	RCP	8	98 15.89 22.07	96 16	16.28	81 15.99 14.81
	MOV	32 118 2.31	5.34	113 2.46	4.21	101 2.46	3.79

This work is sponsored by the French ANR (BioWic project) as well as generous hardware donations from nVidia corporation.


Cyril Bozonnet
email: cyril.bozonnet@univ-grenoble-alpes.fr

Olivier Desjardins

Guillaume Balarac

turbulent, single-or multi-phase flows, and surface wave simulations

Keywords: Outflow, non-reflective boundary, backflow instability

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Due to the finite nature of numerical simulations, it is often necessary to truncate computational domains. This requires imposing artificial boundaries along with the associated mathematical conditions that close the system of equations to be solved. The primary goal of such boundaries is to restrict the computation to a given region of interest without perturbing the solution inside the domain, thereby limiting cost. In the case of outflow boundaries, the flow should be allowed to leave the computational domain in the most natural way possible without undergoing any perturbations that could propagate upstream and thus pollute the upstream solution. Moreover, complex dynamics may occur at the artificial boundary and the flow may contain regions of both outflow and backflow, i.e., regions of flow reversal where the outlet boundary acts as an inlet, potentially polluting the solution [START_REF] Blayo | Revisiting open boundary conditions from the point of view of characteristic variables[END_REF].

The definition of an ideal open boundary condition (OBC) for incompressible fluid dynamic simulations is still an unresolved topic, as demonstrated by Sani and Gresho after the "Open boundary condition minisymposium" [START_REF] Sani | Resumé and remarks on the open boundary condition minisymposium[END_REF],

or by other authors in recent reviews [START_REF] Blayo | Revisiting open boundary conditions from the point of view of characteristic variables[END_REF][START_REF] Bertoglio | Benchmark problems for numerical treatment of backflow at open boundaries[END_REF]. However, one can describe the effect of a non-ideal OBC on a simulation result. In wave-like simulations, the phenomena of wave reflection can create unrealistic flows, instabilities and prevent the flow from reaching a statistical equilibrium over a long computational time [START_REF] Marchesiello | Open boundary conditions for long-term integration of regional oceanic models[END_REF]. In turbulent flows, the presence of backflow can cause the system to experience an uncontrolled growth in kinetic energy, which has for example been evidenced in biofluids simulations [START_REF] Moghadam | A comparison of outlet boundary treatments for prevention of backflow divergence with relevance to blood flow simulations[END_REF].

More generally, the choice of OBC can severely influence the size of the computational domain due to the difficulty of finding a condition that does not durably affect the upstream flow, the most famous example being the impact of the outflow position on a cylinder drag and lift coefficients [START_REF] Persillon | Physical analysis of the transition to turbulence in the wake of a circular cylinder by three-dimensional Navier-Stokes simulation[END_REF].

Indeed, the incompressibility constraint and the unphysical nature of domain truncations may prevent finding a perfect OBC. However, in this paper, we endeavor to present a novel boundary treatment that reduces the error induced by outlet position on severely truncated domains and is stable to backflow, in addition to satisfying the incompressibility constraint.

In the next section, the main types of outflow treatments are discussed.

The new proposed strategy is then presented. In section 3, the numerical implementations of these boundary conditions are presented in the context of a fractional step method with pressure projection method. Section 4 is devoted to single phase test cases, consisting of the Kovasznay flow for measuring spatial convergence, a time-dependent manufactured solution test for measuring temporal convergence, and a flow past a square and a turbulent plane jet to explore the stability and accuracy of the method. Finally, multiphase test cases are considered in section 5, with the convection of a high density droplet, a turbulent swirling liquid jet, and the transport of surface gravity waves. All of the work presented hereafter is applied to outlet boundary con-ditions where the flow is expected to be mostly leaving the computational domain, but it can be applied as well on lateral and inlet boundaries.

Existing methods and present work

Before giving a short review of existing methods we introduce here some useful notations. The computational domain will be refered to as Ω. This domain is bounded by real and artifical (open) boundaries. The firsts ones are denoted ∂Ω d and the laters ∂Ω o . ∂Ω will refer to both types of boundaries, i.e., ∂Ω = ∂Ω d ∪ ∂Ω o . The vector n is defined as the unit normal to those boundaries, always oriented toward the exterior of the domain.

Existing methods

Apart from classical Dirichlet and Neumann conditions, one of the most widely used boundary conditions is the convective boundary condition,

∂φ ∂t + c ∂φ ∂n = 0. (1)
This equation represents the transport of a quantity φ through a boundary of normal n with a phase speed c, where n is the coordinate in the n direction. This condition, known as the Sommerfeld equation, or the radiation condition, is in fact an exact absorbing condition, i.e., specification of the incoming characteristic to zero, for a 1D wave equation with a constant wave speed [START_REF] Blayo | Revisiting open boundary conditions from the point of view of characteristic variables[END_REF]. The most famous choice of phase speed comes from the work of Orlanski [START_REF] Orlanski | A simple boundary condition for unbounded hyperbolic flows[END_REF]: c is computed locally based on known values of φ in the vicinity of the boundary. This solution has been shown to result in a phase velocity close to white noise [START_REF] Durran | Open boundary conditions: fact and fiction[END_REF]. Despite some improvements of Orlanski's method [START_REF] Raymond | A radiation boundary condition for multidimensional flows[END_REF], it seems that no satisfying method has emerged to obtain an accurate estimation of the phase velocity without a priori knowledge of it [START_REF] Blayo | Revisiting open boundary conditions from the point of view of characteristic variables[END_REF][START_REF] Higdon | Radiation boundary conditions for dispersive waves[END_REF].

A general mathematical approach to obtain exact absorbing boundary conditions has been derived [START_REF] Engquist | Absorbing boundary conditions for numerical simulation of waves[END_REF]. However, to our knowledge, no applications of this method to Navier-Stokes equations have been presented, the closest being recent progress on shallow-water equations [START_REF] Blayo | Revisiting open boundary conditions from the point of view of characteristic variables[END_REF]. It has been applied to a 2D wave equation whose coefficients are then identified using the Navier-Stokes equations [START_REF] Jin | A nonreflecting outlet boundary condition for incompressible unsteady Navier-Stokes calculations[END_REF]. It results in a condition similar to Eq. ( 1) with the phase velocity evaluated as the local speed and the presence of a viscous term on the right hand side. More generally, a whole family of OBCs relies on the method of characteristics [START_REF] Blayo | Revisiting open boundary conditions from the point of view of characteristic variables[END_REF].

On the other hand, another type of boundary condition can be directly derived from the Navier-Stokes equation in its weak form [START_REF] Taylor | A note on the imposition of traction boundary conditions when using the FEM for solving incompressible flow problems[END_REF]: the traction boundary condition. It consists of applying a condition on the normal stress at the artificial boundary,

σ • n = (-pI + µ(∇u + ∇u T )) • n = t, (2)
where σ, p, µ and u are the stress tensor, the pressure, the dynamic viscosity, and the velocity, respectively. t is a traction vector that must be prescribed.

No clear guidelines exist for the choice of this vector. The most widespread choice is t = 0, giving the well-known "traction-free" boundary condition [START_REF] Liu | Open and traction boundary conditions for the incompressible Navier-Stokes equations[END_REF][START_REF] Hasan | On the outflow boundary condition for external incompressible flows: A new approach[END_REF]. The traction t has also been computed locally and iteratively [START_REF] Taylor | A note on the imposition of traction boundary conditions when using the FEM for solving incompressible flow problems[END_REF],

based on previous runs on longer domains [START_REF] Sani | Resumé and remarks on the open boundary condition minisymposium[END_REF], or defined analytically with a

Stokes solution [START_REF] Bruneau | Effective downstream boundary conditions for incompressible Navier-Stokes equations[END_REF].

As stated previously, the presence of backflow at an outlet boundary could lead to an instability due to an uncontrolled growth of kinetic energy. To understand it, the energy balance in the overall computational domain, Ω, can be considered [START_REF] Dong | A robust and accurate outflow boundary condition for incompressible flow simulations 37 on severely-truncated unbounded domains[END_REF][START_REF] Dong | An outflow boundary condition and algorithm for incompressible two-phase flows with phase field approach[END_REF],

∂ ∂t Ω 1 2 ρ|u| 2 = - Ω µ 2 D(u) 2 + Ω (ρg + T σ ) • u + ∂Ω d σ • n - 1 2 ρ|u| 2 n • u + ∂Ωo σ • n - 1 2 ρ|u| 2 n • u, (3)
Where ρ is the density, g is the gravity vector, D(u) is the shear rate tensor Following Eq. ( 3), one possible backflow treatment is to ensure that the last term is zero, preventing backflow from causing an unstable growth of kinetic energy. It leads to the following OBC

σ • n = (-pI + µ(∇u + ∇u T )) • n = ρ 2 f (u)n, (4)
with f (u) chosen so that it cancels the last term of Eq. ( 3) in case of backflow, for example

f (u) = (u • n) 2 if u • n < 0, 0 otherwise. ( 5
)
This condition is similar to the stabilized traction-free condition used for single phase flows [START_REF] Dong | A robust and accurate outflow boundary condition for incompressible flow simulations 37 on severely-truncated unbounded domains[END_REF][START_REF] Dong | A pressure correction scheme for generalized form of energy-stable open boundary conditions for incompressible flows[END_REF] and for multiphase flows [START_REF] Dong | An outflow boundary condition and algorithm for incompressible two-phase flows with phase field approach[END_REF][START_REF] Dong | A rotational pressure-correction scheme for incompressible two-phase flows with open boundaries[END_REF]. In case of backflow, the normal stress will compensate the normal influx of kinetic energy, whereas it will vanish in case of outflow. Different forms for f (u) along with other types of backflow treatments have been reviewed [START_REF] Bertoglio | Benchmark problems for numerical treatment of backflow at open boundaries[END_REF].

The traction boundary condition, when used as the stabilized tractionfree condition as in Eq. ( 4), requires the flow to be well-developped before reaching the boundary [START_REF] Taylor | A note on the imposition of traction boundary conditions when using the FEM for solving incompressible flow problems[END_REF]. Several methods have already been proposed to combine stability and accuracy even at high Reynolds number, such as the "convective-like" traction boundary condition [START_REF] Dong | A convective-like energy-stable open boundary condition for simulations of incompressible flows[END_REF],

σ • n = -µD 0 ∂u ∂t + ρ 2 Θ(u • n)((u • n)u + |u| 2 n), (6)
where D 0 is computed using a characteristic velocity, and the function Θ(x) is essentially equal to 1 for negative value of x and 0 otherwise, see [START_REF] Dong | A convective-like energy-stable open boundary condition for simulations of incompressible flows[END_REF] for more details. The value of D 0 is found to have little effect on the overall flow, except on the flow patterns near the outlet boundary. An earlier method developped by Bruneau and Fabrie [START_REF] Bruneau | Effective downstream boundary conditions for incompressible Navier-Stokes equations[END_REF] combines a stabilization to backflow and a non-zero traction,

σ • n = σ ref • n + ρ 2 (u • n) -(u -u ref ), (7)
where the reference values are computed using an analytical solution, or evaluated from known values inside the domain [START_REF] Bruneau | Boundary conditions on artificial frontiers for incompressible and compressible navier-stokes equations[END_REF], and

(u•n) -= max(0, -u•n).
Note that this condition leads to a well-posed problem [START_REF] Bruneau | New efficient boundary conditions for incompressible navier-stokes equations: a well-posedness result[END_REF]. It has, to the best of our knowledge, not been applied to projection methods.

Another potential backflow treatment is to simply force all velocities such that u • n < 0 to zero, thus preventing any influx of kinetic energy. This solution provides energy stability of the system, but we will show in section 5.1 that it can lead to severe inaccuracies in multiphase flows. In [START_REF] Marchesiello | Open boundary conditions for long-term integration of regional oceanic models[END_REF], when the phase velocity is computed as u • n < 0, the use of external data allows to limit the occurence of the backflow instability.

As said in the introduction, the main difficulties encountered by outflow treatments are associated with the proper transmission of perturbations through the artificial boundary, and to the presence of inflow/backflow regions on it. One common way to overcome those issues is to try to dissipate, or damp, the fluctuating energy of the flow before the outlet using artificial zones called sponge layers or nudging layers. Sponge layers consist in the introduction of a dissipative source term in Navier-Stokes equations that becomes stronger when getting closer to the boundary [START_REF] Clément | Coupling of two absorbing boundary conditions for 2D timedomain simulations of free surface gravity waves[END_REF]. Nudging layers consist of the relaxation of the flow towards prescribed external data [START_REF] Marchesiello | Open boundary conditions for long-term integration of regional oceanic models[END_REF].

These solutions are intentionally excluded from our study to focus on the improvements of an accurate OBC.

Finally, most efforts to get non-reflective and accurate boundaries have been focused on convective-like OBCs, often requiring the use of external data that is consistent with the backflow treatment [START_REF] Marchesiello | Open boundary conditions for long-term integration of regional oceanic models[END_REF], whereas traction boundary conditions present an easier way to deal with backflow without the need for external data. As said previously, the stabilized traction-free condition, Eq. ( 4), requires the flow to be well-developped before reaching the boundary [START_REF] Taylor | A note on the imposition of traction boundary conditions when using the FEM for solving incompressible flow problems[END_REF]. Traction boundary conditions have, to our knowledge, never been applied to problems of wave reflections.

Generalized traction boundary condition

We propose a new traction boundary condition, inspired from the Bruneau and Fabrie condition Eq. [START_REF] Orlanski | A simple boundary condition for unbounded hyperbolic flows[END_REF], that combines the two following characteristics.

Firstly, the flow will not be required to be well-developped at the boundary, which will be achieved by applying a non-zero traction at the boundary. Secondly, this OBC will be stable to influxes of kinetic energy due to backflow, which will be achieved by the inclusion of a stabilization term.

We express the traction at the boundary as

(-pI + µ(∇u + ∇u T )) • n = t stab + αt est . ( 8
)
t stab is a numerical treatment to ensure stability to backflow. t est is an estimation of the traction at the outlet boundary and α = [0; 1] is an adjustable parameter. The accuracy of the present boundary treatment will depend on the choice of the last two terms.

Following Eq. ( 3), the stabilization term is taken such that it cancels the term responsible for the backflow instability in case of backflow,

t stab = ρ 2 f (u)n, (9)
with f (u) defined as in Eq. ( 5). Thus, the kinetic energy variation at the open boundary is not equal to zero, as with Eq. ( 4), but depends on the value of αt est . The results presented in this paper show that this novel boundary condition is sufficient to ensure the stability of the system in the presence of backflow at the open boundary. If α is equal to zero, one can see that we recover the stabilized traction-free condition presented in [START_REF] Dong | A robust and accurate outflow boundary condition for incompressible flow simulations 37 on severely-truncated unbounded domains[END_REF].

To obtain the best possible traction estimate we introduce here the general idea behind our work. We propose t est , the estimated traction at the boundary, to be considered as a Lagrangian quantity. Its value can therefore be evaluated using an advection equation,

∂t est ∂t + u ad • ∇t est = 0, (10)
where u ad is an advection velocity that can be computed using an analytical expression, an averaged or a local velocity.

Scope of this work

The previous method to estimate the traction is very general and studying all possible ways to resolve it is beyond the scope of the present paper.

Thus, we restrict our study to a few particular cases. We first assume a one dimensional advection velocity of the estimated traction in the outlet boundary normal direction. Then, we assume a first order explicit temporal resolution of Eq. ( 10) on a cartesian grid. The choice of an explicit resolution is a consequence of the algorithm used to solve the coupling between velocity and pressure, as we detail in the next section. Finally, we use a first order upwind discretization of the spatial term in order to use values inside the computational domain.

The traction estimation is therefore expressed as

t est = φσ BC-1 • n + (1 -φ) σ BC • n , (11)
where the notations BC -1 and BC refer to the point just before the boundary and the boundary point, respectively. φ is an interpolation coefficient computed using numerical parameters and the one dimensional advection velocity. φ can be considered as a CFL condition and therefore has here to be kept in the range [0; 1] as the advection is only done between the boundary point and its closest neighbour.

The previous choices of resolution for Eq. ( 10) are not suitable in case of discontinuities in the traction field. This latter point is limiting in case of multiphase flows due to the effect of surface tension. The presence of a pressure jump can thus deteriorate the traction estimation and create unphysical velocities, or even stability issues. Therefore, in case of multiphase flows we limit our study to high Weber number. A way to get around that difficulty would be to use, for example, a semi-Lagrangian advection method [START_REF] Owkes | A computational framework for conservative, three-dimensional, unsplit, geometric transport with application to the volume-of-fluid (VOF) method[END_REF] for t est . Other aspects may have to be considered, such as the curvature computation in the vicinity of the open boundary, or the density boundary condition. Note that the use of multiphase traction-free condition in phase field method provides a natural way to get around that difficulty as a surface tension term appears in the outlet boundary energy flux [START_REF] Dong | An outflow boundary condition and algorithm for incompressible two-phase flows with phase field approach[END_REF].

For α = 0, the generalized traction boundary condition, Eq. ( 8), reduces to

(-pI + µ(∇u + ∇u T )) • n = ρ 2 f (u)n, (12)
which will be refered to as the stabilized traction-free condition (TF) in the following. TF is the same condition as used in [START_REF] Dong | A robust and accurate outflow boundary condition for incompressible flow simulations 37 on severely-truncated unbounded domains[END_REF]. For α = 1 and φ = 1, Eq. ( 8) reduces to (-pI + µ(∇u

+ ∇u T )) • n = ρ 2 f (u)n + σ BC-1 • n, (13)
which will be refered to as the estimated traction boundary condition (ET).

This condition ressembles the Bruneau and Fabrie condition, Eq. ( 7). The choice φ = 1 raises the question of the dependence of the accuracy to numerical parameters, as the traction at the point just before the boundary may not always be a good estimation. In the final part of the article, we will consider the φ = 1 case, where Eq. ( 8) reduces to

(-pI + µ(∇u + ∇u T )) • n = ρ 2 f (u)n + φσ BC-1 • n + (1 -φ) σ BC • n , (14)
which will be refered to as the convected traction boundary condition (CT).

Note that in the previous three boundary conditions f (u) is computed using Eq. ( 5).

In the rest of the paper we also use classic OBCs, such as the Neumann boundary condition (NM),

∂u ∂n = 0, (15)
or the convective boundary condition (CV),

∂u ∂t + c ∂u ∂n = 0. ( 16
)
As mentioned previously the performance of such condition will be strongly linked to the choice of the convective velocity, which will be detailed later.

Finally, the main objectives of the present paper are, for each of the OBCs under consideration, to give a detailed description of the algorithm allowing to their use in the context of projection methods and VOF/Level Set methods, to demonstrate the importance of bakflow stabilization in single-or multiphase flows, and to demonstrate the stability and accuracy of the nonzero traction methods, such as ET or CT. CT will only be used in the end of the paper, where the level of accuracy obtained with ET is not satisfactory.

Mathematical formulation and algorithms

General framework

Fluid dynamics are governed by conservation laws, forming the Navier-Stokes system of equations. Conservation of mass is, providing that the flow

is incompressible (∇ • u = 0), ∂ρ ∂t + ∇ • (ρu) = ∂ρ ∂t + u • ∇ρ = 0. ( 17
)
Multiphase momentum conservation is written in the framework of the one fluid formulation [START_REF] Tryggvason | Direct numerical simulations of gas-liquid multiphase flows[END_REF]: a single equation with space varying material properties is used to describe the dynamics of both phases. The effect of surface tension is added through a singular force, T σ , acting on the interface,

∂ρu ∂t + ∇ • (ρuu) = -∇p + ∇ • µ ∇u + ∇u T + T σ + ρg. (18)
In absence of phase change, the application of the momentum equation on the interface results in the classical jump condition for normal stress,

[p] Γ = σκ + 2 [µ] Γ n T Γ • ∇u • n Γ , (19)
where σ, κ and n Γ are the surface tension, the curvature and the normal vector to the interface Γ respectively. The notation [.] Γ represents the interfacial jump from liquid to gas.

These equations are solved using NGA, a finite volume, staggered-grid, second order flow solver [START_REF] Desjardins | High order conservative finite difference scheme for variable density low Mach number turbulent flows[END_REF]. Mass conservation, Eq. ( 17), is ensured through an unsplit semi-Lagrangian VOF advection method [START_REF] Owkes | A computational framework for conservative, three-dimensional, unsplit, geometric transport with application to the volume-of-fluid (VOF) method[END_REF]. Momentum conservation, Eq. ( 18), is computed in a way consistent with mass transport and with the presence of interfacial discontinuities [START_REF] Palmore | A volume of fluid framework for interfaceresolved simulations of vaporizing liquid-gas flows[END_REF]. It is to be noted that all OBC methods presented are usable with any sharp interface-capturing method (VOF/Level-Set). For an application to diffuse interface methods, we refer the reader to a work for phase field method [START_REF] Dong | An outflow boundary condition and algorithm for incompressible two-phase flows with phase field approach[END_REF]. Interface boundary conditions, Eq. ( 19), are included in the pressure through the use of the ghost fluid method [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF], with a curvature computed using a least-squares curve fitting method [START_REF] Marchandise | A stabilized finite element method using a discontinuous level set approach for the computation of bubble dynamics[END_REF]. The coupling between velocity and pressure, required due to the incompressibility constraint, is enforced using an incremental pressure projection method [START_REF] Goda | A multistep technique with implicit difference schemes for calculating two or three-dimensional cavity flows[END_REF]. In the following equations, superscripts n and n + 1 refer to previous and new time steps, respectively, whereas subscript k refers to the subiterations of the iterative Crank-Nicolson time advancement scheme [START_REF] Teukolsky | Stability of the iterated Crank-Nicholson method in numerical relativity[END_REF] used in the present solver. Note that we use an implicit resolution of the linearized problem at each subiteration, see [START_REF] Desjardins | High order conservative finite difference scheme for variable density low Mach number turbulent flows[END_REF] for more details. Second order centered schemes are used for spatial discretization on all terms but the convective term at the interface, where a consistent mass and momentum advection strategy is employed (see [START_REF] Palmore | A volume of fluid framework for interfaceresolved simulations of vaporizing liquid-gas flows[END_REF] for more details). In case of single phase flows the same solver is used but the physical properties are taken as constant and surface tension effects are not present. First, a non-solenoidal velocity field, u * k+1 , is computed as

ρ n+1 k+1 u * k+1 -ρ n u n ∆t = -∇p n+1 k -∇ • ρ n u n+1/2 k u n + u * k+1 2 + ∇ • µ n+1 ∇ u n + u * k+1 2 + ∇ u n + u * k+1 2 T + ρ n+1 k+1 g, (20)
where the intermediate velocity field is

u n+1/2 k = 1 2 u n+1 k + u n . (21)
Then, a Poisson equation is solved for the pressure increment Φ n+1 ,

∇ • ∆t ρ n+1 k+1 ∇ Φ n+1 = ∇ • u * k+1 . (22)
Finally, the velocity and the pressure at the next time step are obtained using

Φ n+1 , u n+1 k+1 = u * k+1 - ∆t ρ n+1 k+1 ∇ Φ n+1 , (23)
p n+1 k+1 = p n+1 k + Φ n+1 . (24)
Eqs. ( 20), ( 22) and ( 23)-( 24) are refered to as estimation, projection, and correction. ∆t is the time step size. In the case of multiphase flows, ρ n and µ n+1 values are computed from the old and new VOF field, respectively, whereas ρ n+1 k+1 is computed in a way that ensures consistency between mass and momentum transport [START_REF] Palmore | A volume of fluid framework for interfaceresolved simulations of vaporizing liquid-gas flows[END_REF]. More details can be found in [START_REF] Desjardins | Direct numerical and large-eddy simulation of primary atomization in complex geometries[END_REF].

At all of these steps, boundary conditions have to be provided: velocity boundary condition after estimation and correction, and pressure boundary condition during projection. At the inflow and on the walls, those steps are straightforward and well documented [START_REF] Gresho | Incompressible fluid dynamics: some fundamental formulation issues[END_REF]. For the velocity it simply consists of setting the corresponding values in the velocity vector. As these values will not change during estimation and correction, this step is only necessary after estimation, Eq. ( 20),

u * k+1 ∂Ω d = u n+1 D , (25)
where u n+1 D is an imposed velocity value given by the physics, i.e., inflow or walls. The definition of the pressure boundary condition is directly obtained from the application of Eq. ( 23) on those boundaries,

∂Φ n+1 ∂n ∂Ω d = 0. ( 26
)
The expression of outlet boundary conditions for velocities and pressure at each step of the projection algorithm, resulting in the application of the OBCs presented in section 2, is detailed in the next subsections.

Implementation of convective/Neumann OBC

Neumann (NM) and convective (CV) boundary conditions can be directly used to compute outlet velocities at the estimation step. In the CV boundary condition, Eq. ( 16), a wave velocity c has to be prescribed. In the present work, it is going to be taken as the maximal velocity in the plane just before the exit,

c = c max = max(u * k+1 • n) BC-1 , (27)
or as a theoretical wave speed, if available,

c = c th . ( 28
)
The theoretical expression for the phase velocity will be detailed in the results when used. Except if otherwise stated, the phase velocity will be taken as c = c max . What is of interest here is the definition of the pressure boundary condition that will allow to obtain a solution to the Poisson equation, Eq. ( 22). This comes from the integration of Eq. ( 22) over the computational domain:

∂Ω ∆t ρ n+1 ∇Φ n+1 • ndS = ∂Ω u * k+1 • ndS (29)
Applying Eq. ( 26) will directly lead to the following pressure outlet boundary condition,

∂Ω ∆t ρ n+1 ∂Φ n+1 ∂n ∂Ωo dS = Q in -Q out , (30)
where Q in and Q out are the inlet and outlet flow rates, respectively. Thus, if inlet and outlet flow rates are forced to be the same (including the clipping of negative velocities, as explained in section 2) when considering the application of the velocity OBC and the resolution of the Poisson equation, the pressure outlet boundary condition can simply be a Neumann BC,

∂Φ n+1 ∂n ∂Ωo = 0, (31)
thus ensuring that the integral on the left hand side of Eq. ( 30) is equal to zero. Finally, as the gradient of pressure on all boundaries is equal to zero, there is no need to correct outlet velocities during the correction step. The overall algorithm is presented in algorithm 1.

Input: u n , p n , ρ n in Ω and on ∂Ω 1 Solve Eq. ( 17) using VOF advection → κ n+1 , µ n+1 in Ω 2 for k = 0 to k max -1 do Set all velocities such that u * k+1 • n < 0 to zero in the outlet section 7

Correct outlet flow rate 8 Solve Eq. ( 22) with Eq. ( 26) and Eq. ( 31) → Φ n+1 in Ω 9

Correct velocities Eq. ( 23) and pressure Eq. ( 24) → u n+1 k+1 and p n+1 k+1 in Ω 10 end Output: u n+1 , p n+1 , ρ n+1 in Ω and on ∂Ω Algorithm 1: Algorithm for Neumann and convective OBCs

Implementation of traction-free and estimated traction OBC

The implementation of traction boundary conditions in pressure projection methods has been the subject of many publications in recent years [START_REF] Poux | Improvements on open and traction boundary conditions for Navier-Stokes time-splitting methods[END_REF][START_REF] Liu | Open and traction boundary conditions for the incompressible Navier-Stokes equations[END_REF][START_REF] Dong | A pressure correction scheme for generalized form of energy-stable open boundary conditions for incompressible flows[END_REF][START_REF] Guermond | Error analysis of pressure-correction schemes for the time-dependent Stokes equations with open boundary conditions[END_REF]. See also [START_REF] Angot | Vector penalty-projection method for incompressible fluid flows with open boundary conditions[END_REF] in the context of vector penalty method and [START_REF] Bänsch | A finite element pressure correction scheme for the Navier-Stokes equations with traction boundary condition[END_REF] for an extension to curved artificial boundaries. Note furthermore that all algorithms presented herein may be adapted to velocity correction methods starting from the work presented in [START_REF] Poux | Open boundary conditions for the velocity-correction scheme of the Navier-Stokes equations[END_REF].

The main difficulty is to ensure the validity of the relation,

(-p n+1 k+1 I + µ n+1 (∇u n+1 k+1 + ∇u n+1 k+1 T )) • n = t n+1 (32)
along with the incompressibility constraint at the end of the correction step.

The form of the vector t n+1 will depend on the type of traction boundary condition (TF, ET, or CT, see section 2) and will be explained below. A first strategy consists in simply setting the pressure increment to zero at the outlet [START_REF] Liu | Open and traction boundary conditions for the incompressible Navier-Stokes equations[END_REF], but this strategy is known to limit the order of convergence of the overall method [START_REF] Guermond | Error analysis of pressure-correction schemes for the time-dependent Stokes equations with open boundary conditions[END_REF]. An improvement is found whith an update of the outlet pressure through a rotational pressure correction method [START_REF] Dong | A robust and accurate outflow boundary condition for incompressible flow simulations 37 on severely-truncated unbounded domains[END_REF][START_REF] Dong | An outflow boundary condition and algorithm for incompressible two-phase flows with phase field approach[END_REF][START_REF] Dong | A pressure correction scheme for generalized form of energy-stable open boundary conditions for incompressible flows[END_REF].

It is known with those methods that the use of a rotational pressure correction will significantly improve the convergence order of the overall algorithm.

However, in multiphase flows, there is, to our knowledge, only one example of a rotational pressure-correction [START_REF] Dong | A rotational pressure-correction scheme for incompressible two-phase flows with open boundaries[END_REF], which involves the resolution of a second linear system due to the absence of an analytical solution for the pressure increment. Thus, we choose to employ the method presented in [START_REF] Poux | Improvements on open and traction boundary conditions for Navier-Stokes time-splitting methods[END_REF],

extended for multiphase flow and non-zero traction.

The method presented in [START_REF] Poux | Improvements on open and traction boundary conditions for Navier-Stokes time-splitting methods[END_REF] starts by applying the normal, i.e., perpendicular, projection of the traction boundary condition just after the estimation step, Eq. ( 20), using the available variables, namely p n+1 k and u * k+1 . To simplify the understanding of the method, the coordinate system is taken to be cartesian (x, y, z) with an artificial boundary oriented along x. This first step is then,

-p n+1 k + 2µ n+1 ∂u * k+1 ∂x = t n+1 x , (33)
where t n+1

x is assumed to be known. The final step will be the application of the traction outlet boundary condition on the new variables u n+1 k+1 and p n+1 k+1 ,

-p n+1 k+1 + 2µ n+1 ∂u n+1 k+1 ∂x = t n+1 x . (34)
We are then looking for the pressure increment that will ensure the validity of Eq. ( 34) as well as satisfy the incompressibility condition. We first take the incompressibility condition in the cell just before the boundary,

∇ • u n+1 k+1 = ∂u n+1 k+1 ∂x + ∂v n+1 k+1 ∂y + ∂w n+1 k+1 ∂z = 0, (35)
where we express ∂u n+1 k+1 /∂x using Eq. ( 33) and Eq. ( 34), and ∂v n+1 k+1 /∂y and ∂w n+1 k+1 /∂z using Eq. [START_REF] Bruneau | New efficient boundary conditions for incompressible navier-stokes equations: a well-posedness result[END_REF]. It leads to

∇•u n+1 k+1 = p n+1 k+1 -p n+1 k 2µ n+1 +∇•u * k+1 - ∂ ∂y ∆t ρ n+1 k+1 ∂ ∂y Φ n+1 - ∂ ∂z ∆t ρ n+1 k+1 ∂ ∂z Φ n+1 .
(36) Finally, as the flow is incompressible, the pressure boundary condition is ∂ ∂y

1 ρ n+1 k+1 ∂ ∂y Φ n+1 + ∂ ∂z 1 ρ n+1 k+1 ∂ ∂z Φ n+1 - 1 2µ n+1 ∆t Φ n+1 = ∇ • u * k+1 ∆t , (37)
which is the pressure boundary condition derived in [START_REF] Poux | Improvements on open and traction boundary conditions for Navier-Stokes time-splitting methods[END_REF] adapted to a variable density flow. Previous equations are forming the algorithm used to compute and couple the x-velocity and the pressure at the outflow, and to satisfy exactly the relation Eq. ( 34) along with the incompressibility constraint.

Concerning the tangential components of the velocity, the outflow condition is simply a Neumann condition,

∂v * ,n+1 k+1 ∂x = ∂w * ,n+1 k+1 ∂x = 0. ( 38
)
This choice, rather than the use of a constraint on the tangential traction value, is motivated by the well-known fact that a tangential traction-free condition is not compatible with a parallel flow [START_REF] Leone | Finite element simulations of steady, twodimensional, viscous incompressible flow over a step[END_REF] and by the fact that several results are reported as better with a Neumann condition on tangential velocities rather than a tangential traction condition, even with non-zero traction [START_REF] Gresho | Incompressible fluid dynamics: some fundamental formulation issues[END_REF].

It should be noted that the pressure boundary condition, Eq. ( 37), is only

valid if t n+1
x does not change between the estimation and correction steps.

Otherwise, any change will have to be taken into account into the pressure OBC, Eq. ( 37). Thus, the traction t n+1

x can be given depending on the type of open boundary condition. For TF, it is

t n+1 x = ρ n 2 f (u n ), (39)
and for non-zero traction conditions (ET and CT),

t n+1 x = ρ n 2 f (u n ) + t est,n+1 x . ( 40
)
The density is taken at the previous time step to be coherent with the choice of the velocity. The backflow stabilization is thus not instantaneous but delayed by one time step. As stated previously, the estimated normal traction is computed using interior values and at the previous iteration to ensure the validity of the pressure boundary condition. For ET,

t est,n+1 x = -p + 2µ ∂u ∂x n BC-1 , (41)
and for CT,

t est,n+1 x = φ -p + 2µ ∂u ∂x n BC-1 + (1 -φ) -p + 2µ ∂u ∂x n BC , (42)
with φ to be prescribed later. The overall algorithm is presented in algorithm 2.

Input: u n , p n , ρ n in Ω and on ∂Ω 1 Solve Eq. ( 17) using VOF advection → κ n+1 , µ n+1 in Ω 2 for k = 0 to k max -1 do Apply Eqs. ( 25)-( 33) and ( 38) on u * k+1 6

Solve Eq. ( 22) with Eq. ( 26) and Eq. ( 37) → Φ n+1 in Ω 7 Correct velocities Eq. ( 23) and pressure Eq. ( 24

Kovasznay flow

The Kovasznay flow is a steady state flow used to mimic the flow behind a cylinder [START_REF] Kovasznay | Laminar flow behind a two-dimensional grid[END_REF]. This configuration is a 2D domain, periodic along the vertical axis, with an inflow on its left boundary and an OBC on its right. The analytical solution of the Kovasznay flow is given by Hereafter we study the effect of the OBC choice on the error compared to the theoretical solution. Considered OBCs are NM, TF and ET. CV is intentionally excluded from this test case to avoid any discussion on the choice of the convective velocity at this point.

u = 1 -e λx cos(2πy), ( 43
) v = λ 2π e λx sin(2πy), (44)
p = 1 2 1 -e 2λx , (45)
In a first comparative test, the domain length is kept constant with L x = 4.5 and the mesh is progressively refined in order to check the convergence of the error depending on the type of OBC. In figure 2a) we show the evolution of the L 2 error norm of the x-velocity and the pressure, for differents OBCs and depending on mesh resolution. One can first observe two differents behaviors:

for coarser meshes, the same level of error is obtained for all three OBCs, which decreases with mesh resolution (with order 2, i.e., the order of used numerical methods). For finer meshes, and for NM and TF OBCs, the error progressively saturates at a constant value, indicating that outflow error is dominating. Note that this deviation occurs later for TF than for NM and stabilizes also at a lower value, meaning that TF gives a lower error than NM on that test case. On the other hand, with ET, no deviation is observed from the second order slope, meaning that in that range of mesh resolutions, the error due to the outflow is never dominating. With finer meshes and ET, one will necessarily observe a saturation of the error as the choice of the estimated traction is not perfect. Note that one can also compute the estimated traction using the analytical solution [START_REF] Dong | A robust and accurate outflow boundary condition for incompressible flow simulations 37 on severely-truncated unbounded domains[END_REF], which is not possible in real flow simulations. Note that second order convergence is also obtained in

L inf error norm.
In a second comparative test, the mesh is kept constant (∆ = 1/80) and the domain is progressively truncated. Similarly to the previous test, we show on figure 2b) the evolution of the L 2 error norm of the x-velocity and the pressure, for differents OBCs and depending on the position of the artificial boundary. It is seen in figure 2b) that on a sufficiently long domain all OBCs produce the same level of error. It is also seen that with NM the truncation of the domain has a much stronger effect than with other OBCs.

The result is, for the range of L x considered here and using NM, barely independent of the artifical boundary position. This point is improved with TF, which provides a better independence of the result with the position of the outlet. With ET the result is independent of L x for a large range of domain size, even though a small increase of L 2 (u) is noticeable. Note that all OBCs are stable for the smallest domain (L x = -0.1), where the outflow boundary is located in a recirculation zone, which is not possible without backflow treatment such as clipping or stabilization. But, in case of the durable presence of a backflow, even the traction condition will not provide a perfectly accurate solution as the stabilization term only makes sense in terms of kinetic energy conservation. In order to examine how the order of spatial convergence is deteriorated whith domain truncation, we show in figure 3, two additional convergence tests with domains shorter than the one used for the test presented in figure 2a). For L x = 0.5, see figure 3a), second order convergence is only obtained for L 2 (u) using ET and for a shorter range of resolution than previously. For L x = -0.1, see figure 3b), the outlet boundary is located in the recirculation zone, and no clear convergence of the error with resolution can be observed, even using ET. For all resolutions and domain lengths tested, ET has the lowest level of error. The same results are obtained on L inf error norm.

10 -2 10 -1 10 -4 10 -3 10 -2 10 -1 ∆ Error norm L 2 (u) L 2 (p) a) 10 -2 10 -1 10 -4 10 -3 10 -2 10 -1 ∆ Error norm L 2 (u) L 2 (p) b)
Those tests demonstrate the interest of the non-zero traction OBC on a steady state problem in terms of error level and independence to outlet position.

Time-dependent manufactured solution

In order to study the temporal order of convergence of the proposed method, we use the time-dependent manufactured solution of [START_REF] Dong | A robust and accurate outflow boundary condition for incompressible flow simulations 37 on severely-truncated unbounded domains[END_REF], u = 2 cos(πy) sin(πx) sin(t), [START_REF] Evrard | A multi-scale approach to simulate atomization processes[END_REF] v = -2 cos(πx) sin(πy) sin(t),

p = 2 sin(πx) sin(πy) cos(t), [START_REF] Munk | The solitary wave theory and its application to surf problems[END_REF] which satisfies the incompressibility condition (∇ • u = 0). In order to satisfy Eq. ( 18), unsteady body forces have to be added to the Navier-Stokes equations.

The computational domain is two-dimensional, of size 0 x 2 and -1 y 1, with 256 uniform cells in both directions. Eqs. ( 46) and [START_REF] Bruneau | Far field boundary conditions for incompressible flows computation[END_REF] are enforced as Dirichlet boundary conditions on three boundaries of the computational domain, whereas the traction condition Eq. ( 32) is used on the last one. Similarly to [START_REF] Dong | A robust and accurate outflow boundary condition for incompressible flow simulations 37 on severely-truncated unbounded domains[END_REF][START_REF] Dong | A rotational pressure-correction scheme for incompressible two-phase flows with open boundaries[END_REF][START_REF] Poux | Improvements on open and traction boundary conditions for Navier-Stokes time-splitting methods[END_REF], the right hand side of Eq. ( 32) is computed The results are shown on figure 4. One can see that the error norm convergence for all flow variables is approximately of second order until a progressive saturation of the temporal error by the spatial error. Note that the convergence of the error seems to be faster for the velocities than for the pressure, which may be due to the presence of splitting errors [START_REF] Poux | Improvements on open and traction boundary conditions for Navier-Stokes time-splitting methods[END_REF]. A rotational pressure correction would be a solution to resolve this discrepancy [START_REF] Poux | Improvements on open and traction boundary conditions for Navier-Stokes time-splitting methods[END_REF][START_REF] Dong | A robust and accurate outflow boundary condition for incompressible flow simulations 37 on severely-truncated unbounded domains[END_REF], but for the reason cited in section 3.3 we chose not to use this strategy. Note that the same results are obtained in L inf error norm. In agreement with [START_REF] Poux | Improvements on open and traction boundary conditions for Navier-Stokes time-splitting methods[END_REF], these results suggest that the present algorithm for the implementation of traction conditions does not deteriorate the order of temporal convergence.

The choice of manufactured solution and computational domain we made involves u • n = 0 and p = 0 at the boundaries (same test case as [START_REF] Dong | A robust and accurate outflow boundary condition for incompressible flow simulations 37 on severely-truncated unbounded domains[END_REF]). In order to prove that the results are not affected by these choices, we use the same manufactured solution on a computational domain shifted along x, i.e., The results are shown in figure 5a). Similarly to the previous results, second order temporal convergence is obtained, though with an error slightly higher than for the previous computational domain. This demonstrates once again that the present algorithm for traction condition does not deteriorate the order of temporal convergence. When we add the stabilization term to the analytical traction value, as done for the results presented in figure 5b), we observe degraded temporal convergence leading to an error plateau. We emphasize that this plateau is fully expected, as the stabilization term is purely ad-hoc and represents a numerical error when added to the analytical traction at the open boundary.

L 2 (u) L 2 (v) L 2 (p)

Flow around a square

We now compare the different OBCs on an unsteady case: the flow over a two-dimensional square. This test case presents two main interests from the point of view of OBC performance. Firstly, we study their ability to convect the vortices generated by the von Kármán instability through the 21 We observed that if no backflow stabilization is taken into account, vortices were gathering on the artificial boundary, finally leading to the blow-up of the simulation due to the backflow instability. Note that, although the vortices exit seems unnatural, once the stabilization term is included the code remains perfectly stable to backflow at the outlet boundary. On the other hand, with ET, no vortex sticking is observed and vortices simply cross the boundary with barely any deformation.

10 -4 10 -3 10 -2 10 -4 10 -3 10 -2 10 -1 ∆t Error norm L 2 (u) L 2 (p) L 2 (v) a) 10 -4 10 -3 10 -2 10 -4 10 -3 10 -2 10 -1 ∆t Error norm L 2 (u) L 2 (p) L 2 (v) b)
We now propose a more quantitative comparison between different OBCs through the study of the aerodynamic quantities. To avoid any confinement effect and any impact of the inflow position, the domain is this time of size Figure 7 presents the evolution of differents aerodynamic quantities as a function of the position and the type of open boundary. We also included the results of two recent publications using the stabilized traction-free condition [START_REF] Dong | A robust and accurate outflow boundary condition for incompressible flow simulations 37 on severely-truncated unbounded domains[END_REF] and the traction-free condition [START_REF] Poux | Improvements on open and traction boundary conditions for Navier-Stokes time-splitting methods[END_REF]. In order to simplify comparison with other publications, the evolution of the aerodynamic quantities is presented in terms of error relative to the value obtained on the longest domain. Figure 7a) shows the evolution of the mean drag coefficient, figure 7b) shows the evolution of the r.m.s lift coefficient and figure 7c) shows the evolution of the Strouhal number associated with the vortex shedding frequency. The results obtained with ET keep a correct behavior even on the smaller domains by exhibiting the lowest variation of aerodynamic quantities as a function of the outflow position. This improvement is largely explained by the fact that the underlying assumption of well-developped flow associated with the tractionfree condition, stabilized or not, is no longer required with ET.

Turbulent plane jet

To finally assess the stability and the accuracy of the proposed boundary condition we study the spatial evolution of a turbulent plane jet. The configuration, the expression of the analytical inlet velocity profile and the choice of parameters are the same as in Da Silva & Metais [START_REF] Da Silva | On the influence of coherent structures upon interscale interactions in turbulent plane jets[END_REF] (case refered to as "DNS2" in their original paper). The numerical domain is a 3D domain of size 0 x 12.4h, -6h y 6h and -1.6h z 1.6h, where h is the jet width. The inlet boundary is located at x = 0 and the outflow at x = 12.4h.

The other boundaries are periodic. The domain is discretized with a uniform cell size ∆x = ∆y = ∆z = 0.04h. The constant time step is ∆t = 0.02.

In a first study the Reynolds number based on the jet width, Re = (U i -

U f f )h/ν
, is taken equal to 3000, with U i the jet centerline inlet velocity and U f f the inlet far-field velocity. The isocontours of positive Q-criterion [START_REF] Dubief | On coherent-vortex identification in turbulence[END_REF] are shown on figure 8, using TF (a) and ET (b). On both figures one can see the spatial development of the jet, initiated by the apparition of successive Kelvin-Helmholtz rolls that are then connected by the apparition of vortices in the streamwise direction. When reaching the outlet boundary the flow is fully tridimensional. On figure 8a) one can see the dramatic effect of TF on the exit of the vortices. A part of them sticks to the boundary and is prevented to leave the domain. On the other hand, using ET, no vortex sticking is observed and the vortices are crossing the boundary with barely any deformations, see figure 8b). With CV, the vortices exit looks very similar to the one seen using ET (results not shown). Note that using TF we had to reduce the time step size in order to obtain a stable simulation.

This may be due to the vortex sticking phenomena coupled to the delayed backflow correction of one time step. Considering the poor qualitative result seen on figure 8a) and the need to decrease the time step size, we therefore exclude TF from the following analysis.

The jet exhibits a self-similar behaviour in its downstream region [START_REF] Da Silva | On the influence of coherent structures upon interscale interactions in turbulent plane jets[END_REF]. In this region, several quantities computed from the time averaged velocity field are evolving linearly with the downstream distance. We choose here to use the centerline velocity, U c = u(x, y = 0) , and the jet half-width, δ 1/2 , defined as the y-location where the velocity is equal to half of the centerline velocity, i.e., u(x, y

= δ 1/2 ) -U x,∞ = 0.5 (U c -U x,∞ ), with U x,∞ = u(x, y = ∞) ,
the far-field velocity. Those quantities follow the following relationships [START_REF] Gutmark | The planar turbulent jet[END_REF]:

δ 1/2 h = K u1 x h + K u2 , (49)
and

U i -U f f U c -U x,∞ 2 = C u1 x h + C u2 . (50)
We plot those quantities, along with their linear relations, on figure 8c) and 8d) using ET and CV. Using ET, in addition to provide a natural exit of the vortices as well as the stability of the simulation, the self-similar region is barely disturbed by the presence of the open boundary. With CV both of the self-similar quantities are strongly affected by the presence of the outflow.

Note that the slopes of the linear relations are the same as in Da Silva & Metais [START_REF] Da Silva | On the influence of coherent structures upon interscale interactions in turbulent plane jets[END_REF].

To demonstrate that the proposed boundary treatment is stable for highly turbulent flows, we study the large eddy simulation of the turbulent plane jet at Re = 30000. The sub-grid stresses are estimated with a dynamic Smagorinsky eddy viscosity model using Lagrangian averaging to compute the dynamic coefficient [START_REF] Meneveau | A lagrangian dynamic subgridscale model of turbulence[END_REF]. The computational domain is now larger in the vertical direction to account for the entrainment induced by the jet, i.e., -8h y 8h. The isocontours of positive Q-criterion are shown on

[(U i -U f f )/(U c -U x,∞ ] 2 ET CV M odel d)
t(U i -U f f )/h V ρU 2 dV /(2ρ(U i -U f f ) 2 ) b)

Multiphase test cases

We now turn our attention to multiphase flows. In this section, the importance of backflow stabilization is demonstrated using first a single drop advection test case, then a turbulent swirling jet flow simulation. Finally, we demonstrate the improvements obtained using CT on a problem of surface waves reflection. On the top pictures of figure 10, one can observe that with CV OBC the drop is flattening on the boundary and no liquid is exiting the domain. On the other hand, using ET, the drop is completely going out with minimal deformation. The reason for these completely different behaviors lies in the fact that the incompressibility condition requires the outlet flow rate to be equal to the inlet flow rate, in this case zero. Thus, the only way for the drop to exit is to allow backflow. We see here one strong limitation of the clipping strategy, which severely affects the flow by preventing the drop from going out, though it provides unconditional stability. On the other hand, once the stabilization term is included, ET and TF are perfectly stable to backflow as can be seen in figures 10f-g-h).

Turbulent swirling jet

To show the importance of backflow stabilization in a more realistic case, we present a simulation of turbulent swirling jet. As shown in figure 11, a turbulent liquid jet exits from a nozzle located on the left of the domain.

The jet then develops into a conical shape and becomes subject to different interfacial instabilities leading to its atomization. The outflow is located on the right of the domain (colored in pink), whereas all lateral boundary conditions are periodic. All physical properties, injection parameters and geometric characteristics are the same as in [START_REF] Evrard | A multi-scale approach to simulate atomization processes[END_REF]. The domain is discretized with a 200 × 400 × 400 cartesian grid and the simulation is advanced with a CFL number of 0.8.

On figures 11a) and 11b), the liquid-gas interface colored by the axial velocity is shown after a simulation time tU/D = 11, where U is the bulk injection velocity, and D the external diameter of the injector. On figure 11a), the result with CV is shown, with the phase velocity computed using Eq. ( 27). One can see that since this outflow treatment does not allow backflow, some of the liquid is prevented from going out and "splashes" on the exit plane. Figure 11b) shows the result using the ET OBC, and in this case the liquid is not blocked on the exit plane. This obviously has a large impact on the capability to reach long term simulations of such atomizing liquid jets. With a boundary condition that does not allow backflow, the simulation time is obvisouly limited by the length of the domain, which is not the case with a stabilized traction boundary condition. It should be noted here that CV may be replaced by NM and ET by TF for the same results on the liquid exit. The difference between TF and ET will lie in the speed of the droplets in the vicinity of the outlet and in the behavior of the vortices exiting the domain, as already discussed using the test case of the flow around a square.

Although here traction conditions are used only as exit conditions, it should be noted that they may also be used as lateral boundary condition [START_REF] Dong | A convective-like energy-stable open boundary condition for simulations of incompressible flows[END_REF][START_REF] Bruneau | Far field boundary conditions for incompressible flows computation[END_REF]. One can also imagine replacing the wall used around the liquid injector by an open traction boundary condition in order to get a more realistic representation of such jets by allowing the development of a "natural" gas co-flow. We will investigate that point in future works.

Surface gravity waves

We finally evaluate the ability of the different OBCs to evacuate a surface wave without reflection. As said in the introduction, wave reflection is a problem of critical importance in ocean modeling as it prevents the convergence of flow statistics and may create unrealistic flows [START_REF] Marchesiello | Open boundary conditions for long-term integration of regional oceanic models[END_REF].

The test case is set up using solitary wave theory [START_REF] Munk | The solitary wave theory and its application to surf problems[END_REF]. The interface height is defined as

η(x) = A 0 sech 2 3A 0 4h 0 3 x , (51)
with A 0 the initial height of the wave and h 0 the water depth. The initial velocity is defined as u = (u(x), 0, 0) where In figure 12c), we use the theoretical wave speed as the convective velocity in CV. One can see that the reflection is much lower in amplitude but creates an increase of the mean liquid level, seen also in figure 13. Thus, even with the best choice of the wave speed, a convective condition is not able to evacuate a soliton out of the domain without reflection.

u(x) = η(x) |g| (h 0 + A 0 ) h 0 + η(x) + U in , (52)
We now focus our study on the use of a traction condition. First, it should be noted that the use of TF is impractical for such simulations. The pressure being hydrostatic in the domain, using a traction-free condition will impose a pressure close to zero at the outlet (velocity gradients being small far from the soliton), thus resulting in a strong suction of the flow which rapidly propagates up to the inlet. On figure 12b), we show the result using the ET. One can observe that the reflection is almost suppressed but that small waves propagate upstream. Thus, such an arbitrary choice of the estimated traction (computed at the point just before the boundary at the previous time step) is in fact even better than the best choice of a convective OBC.

However, two points have to be emphasized. Firstly, ET is not perfectly non-reflective. Secondly, we have observed a dependence of its performance to numerical parameters such as the time step or cell sizes.

The reason for the last two points has in fact already been explained in section 2. Considering the estimated traction as a Lagrangian quantity, the present choice is not optimal. Indeed, taking the estimation of the traction at the point before the boundary and at the previous time step is not always a good choice depending on numerical parameters. A more accurate choice can be found considering the convected traction boundary condition (CT).

The theoretical wave speed, c th , is taken as the advection velocity. The interpolation coeficient used in Eq. ( 14) is therefore φ = c th ∆t/∆x. The result using this approach is shown on figures 12d) and 13).

One can see that CT, our new OBC, is now perfectly non-reflective. It must be emphasized that this result is now independent of the time step, the mesh size, and is also independent of the OBC position.

Additional remark

In most of the present paper, we use an estimation of the traction at the point just before the boundary, which is easy to define on a structured mesh. The extension to fully unstructured meshes is possible thanks to the Lagrangian estimation of the traction introduced previously, as done in Eq. ( 14). On such meshes, one has to define an advection velocity and then perform a semi-Lagrangian interpolation of the traction field at the location of interest to obtain the traction estimate. In case of curved boundaries, one has to use differential geometry to complete the pressure boundary condition, Eq. ( 37), as done in [START_REF] Bänsch | A finite element pressure correction scheme for the Navier-Stokes equations with traction boundary condition[END_REF].

Conclusion

We have presented a comparison between several outlet boundary treat-

  and T σ represents surface tension forces. It results that the rate of change of kinetic energy is controlled by viscous dissipation (exchange with internal energy), gravity (exchange with potential energy), surface tension (exchange with surface energy) and by two surface terms. The first one is expressed on ∂Ω d , the Dirichlet boundaries, where variables are known. The second surface term is expressed on ∂Ω o , the outflow boundary, where all variables have to be computed. In case of backflow, the convective part of this term becomes positive and can lead to a global increase of kinetic energy, leading to the instability of the system.

3ρ n+1 k+1 4

 4 ComputeSolve Eq. (20) → u * k+1 in Ω 5Apply Eq. (25) and Neumann or convective OBC on u * k+1 6

3ρ n+1 k+1 4

 4 ComputeSolve Eq. (20) → u * k+1 in Ω 5

) → u n+1 k+1 and p n+1 k+1 in Ω 8 Apply

 8 Eqs.[START_REF] Gresho | Incompressible fluid dynamics: some fundamental formulation issues[END_REF] and[START_REF] Bänsch | A finite element pressure correction scheme for the Navier-Stokes equations with traction boundary condition[END_REF] on u n+1 k+1 9 endOutput: u n+1 , p n+1 , ρ n+1 in Ω and on ∂Ω Algorithm 2: Algorithm for traction OBCs4. Single phase test casesThe improvements obtained using our novel outlet treatment are first illustrated on single phase test cases. The first test case, the Kovasznay flow, allows to see the spatial order of convergence of the overall method, while the second test case, a time-dependent manufactured solution test, allows to study the temporal order of convergence of the present algorithm. The third test case, the flow around a square, shows both qualitative and quantitative improvements thanks to ET. The last case, a turbulent plane jet, shows the stability and accuracy of ET in the presence of a fully turbulent flow.

Figure 1 :

 1 Figure 1: Streamlines of the Kovasznay flow

Figure 2 :

 2 Figure 2: Error levels for different OBCs: a) Mesh convergence, b) Effect of domain truncation. Continuous line: ET -Dashed line: TF -Dotted line: NM

Figure 3 :

 3 Figure 3: Error levels for different OBCs and two domain lengths: a) L x = 0.5, b) L x = -0.1. Continuous line: ET -Dashed line: TF -Dotted line: NM. The orange dashed line in both plots shows second order spatial convergence.

  using the manufactured solution. Imposing the analytical traction value at the open boundary is needed to avoid the domination of spatial error over temporal error. No stabilization term is included as it becomes meaningless when one imposes the exact traction at the open boundary, i.e., when the open boundary is transformed into a Dirichlet boundary. The initial velocity field is set to zero, in agreement with the manufactured solution. For this test case we use ρ = 1 and µ = 0.01. The simulation is advanced in time with a fixed time step, ∆t, to be specified, until a fixed final time t f = 0.5. The L 2 error norm on different flow variables at this final time is computed with respect to the manufactured solution. The test is then repeated with various time steps.

Figure 4 : 1 .5 and -1 y 1 .

 411 Figure 4: Method of manufactured solution. Convergence of the error level with respect to the time step. The dashed black line shows second order convergence.

Figure 5 :Figure 6 :

 56 Figure 5: Method of manufactured solution for the shifted domain. Convergence of the error at the final time with respect to the time step, a) imposing the analytical traction at the open boundary, and b) imposing the analytical traction and the stabilization term at the open boundary. The dotted and dashed black lines show first and second order convergence, respectively.

  L the position of the outlet boundary and H the size of the square located at (0, 0). The Reynolds number Re = ρU H/µ is now equal to 100. The domain is discretized with a uniform cell size ∆ = H/40 in the sub-domain -10H x L and -4H y 4H, to avoid any loss of resolution of the vortices in the wake of the obstacle, and is then progressively stretched up to the top and bottom boundaries with a constant stretching ratio of 1.05. The time step is chosen such that the CFL number stays equal to 1. The aerodynamic forces are directly integrated on the surface of the obstacle.

Figure 7 :

 7 Figure7: Impact of the outflow position and type on the aerodynamic quantities: a) mean drag coefficient, b) r.m.s lift coefficient, c) Strouhal number. On all plots, the vertical axis is the error compared to the value obtained on the longest domain, whereas the horizontal axis is the distance between the square and the outflow position L normalized by the size of the square H. The reference results included are the results from Dong et al[START_REF] Dong | A robust and accurate outflow boundary condition for incompressible flow simulations 37 on severely-truncated unbounded domains[END_REF] and Poux et al[START_REF] Poux | Improvements on open and traction boundary conditions for Navier-Stokes time-splitting methods[END_REF].

Figure 8 :

 8 Figure 8: (Color online) Turbulent plane jet at Re = 3000. Positive Q-criterion isocontours (20 isocontours from Q= 0.25 to Q = 100) at t(U i -U f f )/h =166, using TF (a) and ET (b). c) Evolution of the jet half-width with downstream distance using ET and CV. Model computed using Eq. (49) with K u1 = 0.089 and K u2 = -1. d) Evolution of the centerline jet velocity with downstream distance using ET and CV. Model computed using Eq. (50) with C u1 = 0.165 and C u2 = -2.

Figure 9 :

 9 Figure 9: Large eddy simulation of a turbulent plane jet at Re = 30000. a) Positive Q-criterion isocontours (10 isocontours from Q = 0.25 to Q = 100) at t(U i -U f f )/h = 400 using ET. b) Temporal evolution of the normalized kinetic energy integrated over the domain using ET.

Figure 10 :

 10 Figure 10: Water drop advection in a domain without inlet. a)-d) : CV ; e)-h) : ET. From left to right, all figures are separated by a time interval of 0.6D/U l .

Figure 10

 10 Figure10shows velocity vectors along with the liquid-gas interface during the advection of the drop towards the outlet for two types of boundary conditions. On the top row of images the result are obtain with the CV, and on the bottom row of images the result are obtained with ET. Note that, on this test case, one can replace CV by NM and ET by TF, for the same qualitative result.

Figure 11 :

 11 Figure 11: (Color online) Turbulent swirling jet test case. Liquid-gas interface colored by axial velocity shown at a time tU/D = 11 (dark blue: u = 0, dark red: u = 3U , with U the bulk injection velocity). a) result using CV, b) result using ET.

  and U in is the inflow velocity. The computational domain is two-dimensional, of size -60h 0 x 20h 0 and -h 0 y 4h 0 with symmetry boundary conditions along y, a constant velocity inflow u = U in at x = -60h 0 and the OBC at x = 20h 0 . Air/water conditions are used for the choice of physical properties. This setup results in the transport of a soliton from the position x = 0 to the OBC at a constant phase velocityc th = |g| (h 0 + A 0 ) + U in .For all cases presented below parameters are chosen as A 0 = 0.005, h 0 = 0.01,U in = 0.07. The domain is discretized using a uniform Cartesian mesh with ∆x = ∆y = 5 × 10 -4 . The solution is advanced using a time step size ∆t = 1 × 10 -3 . On figure 12 are presented the space-time plots of the interface height for 4 different OBCs along with, on figure 13, the interface height signals at a position x = 10h 0 . On figure 12a), the result with CV is shown. As in previous tests, the wave speed is taken as c = c max . One can first see a transient phenomenon at the initialization which causes the emission of perturbations towards the left of the domain and the height of the wave to slightly decrease. Since the inflow is located at x = -60h 0 , none of the results presented herein are affected by the reflexion of these initial perturbations on the inflow. The reason for these perturbations is an initial adjustment due to the discrete approximations of the continuous solution [4]. Following this initial transient, the soliton travels towards the OBC at a constant speed c th . Once the wave reaches the artificial boundary, it completely crashes on the boundary and a large part is reflected in the domain in a succession of smaller waves, forming a reflection cone. This result is not a surprise given the inapropriate choice of the convective velocity.

Figure 12 :

 12 Figure 12: (Color online) Transport of a surface gravity waves through an OBC. a) CV with c = c max , b) ET, c) CV with c = c th , d) CT with c = c th . The color indicates the liquid height.

Figure 13 :

 13 Figure 13: Interface height signals at a fixed position x = 10h 0 for 4 different OBCs.

  ments on single and multiphase test cases along with their numerical implementation in the context of fractional step methods. One major difference between these open boundary conditions lies in the backflow treatment. The implementation of backflow clipping associated with Neumann or convective open boundary conditions, while providing unconditional stability, can have a strong effect in the simulation of multiphase flows. On the other hand, stabilized traction conditions are perfectly suited to resolve this issue. The main drawback of the traction-free condition lies in its underlying assumption of well-developed flow that is not suited for severely truncated domains or high Reynolds number flows. To overcome this issue, an open boundary condition combining stabilization to backflow and space and time varying estimated traction is proposed, allowing stable and accurate simulations for turbulent and multiphase flows. This estimated traction is considered as a Lagrangian quantity, which allows to use it as a non-reflective artifical boundary for surface waves simulations. This work shows that traction conditions have the potential to resolve most of issues related to outflow treatment. They might also be used as lateral or inlet boundary conditions and allow a considerable reduction in the cost of numerical simulations, as we will explore in future work. The very general form under which the estimated traction is introduced also opens the way to a study of the effect of different advection methods on the accuracy of traction boundary conditions.

Acknowledgements

We acknowledge Prof. Eric Blayo, for pointing to our attention the work of [16] and for several discussions, as well as Dr. Alexandre Poux and Prof. Mahdi Esmaily for fruitful discussions. This work was sponsored by the Office of Naval Research (ONR) as part of the Multidisciplinary University Research Initiatives (MURI) Program, under grant number N00014-16-1-2617, and by the IDEX UGA "International Strategic Partnerships" program. A part of the simulations were performed using HPC resources from GENCI-CINES (Grant No. A0072A00611).