A moment-of-ﬂuid method for resolving ﬁlamentary structures using a symmetric multi-material approach

Multiphase ﬂows have implications in many areas of engineering. The moment-of-ﬂuid (MOF) method is an interface capturing method using both volume fraction and centroid within a cell for interface reconstruction. A symmetric approach to reconstruct thin structures is presented. Also called ﬁlaments, these subcell characteristics involve multi-material reconstruction. In addition, a new optimisation algorithm is presented using a bisection method without any necessary initial condition. Using a Lagrangian approach for dynamic cases, no restrictions are imposed on timestep. The new method is validated using several benchmark cases that have been studied extensively in the literature. A near quadratic order of convergence and high accuracy is achieved while maintaining an acceptable runtime.  2023 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons .org /licenses /by /4 .0/).


Introduction
Multiphase flow modelling is crucial in several real-life examples such as wave breaking, water splashing or bubbles.It is also important in industrial applications from oil-and-gas transportation to inkjet printing.Its modelling requires an accurate representation of the interface between two or more fluids.In addition, it is challenging to resolve thin filaments during the breakup process.
Several techniques for representing interfaces have been developed over the years [1].These fall into two broad categories: interface tracking and interface capturing.Interface tracking aims to track a set of points representing the interface using the associated velocity field [2][3][4][5].Its robustness and simplicity of implementation makes it easily accessible.However, this may not be the case when there is a large deformation of the interface.
The level set method is an example of an interface capturing technique.This method uses a smooth function to describe the sharp interface [6].The zero level set of the function defines the interface.The accuracy and robustness makes it useful for complex multiphase flows, yet in most engineering problems, the lack of mass conservation of this technique makes it undesirable.Improvements have been made on this issue using a conservative level set method [7,8], which has also been extended recently for non-Newtonian multiphase flows [9].Another interface capturing technique is the volume-of-fluid (VOF) method which was initially developed by Hirt and Nichols [10].The interface, either horizontal or vertical known as Simple Line Interface Calculation (SLIC), is defined by the volume fraction of surrounding grid cells.In a subsequent development, Young's method introduced an improved orientation to the interface known as Piecewise Linear Interface Construction (PLIC) [11].Both VOF methods are subject to natural diffusion and artificial surface tension causing the separation of forming filaments and exhibit large errors.In addition, only one interface is able to be reconstructed in a cell, which means that structures thinner than a cell size cannot be resolved.More sophisticated VOF approaches have been developed using a parabolic reconstruction [12,13] or trying to resolve filaments [14].Methods coupling the level set and VOF methods have also been developed [15].
The latest advancement in the evolution of the VOF method is the moment-of-fluid (MOF) method.The MOF method uses both the volume fraction and its centroid to reconstruct the interface thereby increasing the accuracy in interface orientation.The error in interface reconstruction is improved when compared to a standard VOF method and the MOF method possesses superior mesh convergence properties.The MOF method is able to reconstruct a piecewise linear interface without using information from neighbouring cells.In addition, MOF can be implemented with ease for general polyhedral cells.However, MOF can be time-consuming as the computational bottleneck is an optimisation algorithm which is required in order to reconstruct the interface in each cell.
The initial MOF method was presented in 2005 using an optimisation algorithm [16].While conserving mass in each cell, the best approximation is to find the normal to the interface that minimises the distance between the reference and reconstructed centroid.In a subsequent development, an analytical solution which avoids the need to employ an optimisation algorithm was found by Lemoine et al. [17].However, this is restricted to rectangular cells.Further work has been performed with the MOF method using more complex approaches.Multi-material reconstruction is facilitated compared to VOF methods when three or more materials are present [18,19], under-resolved filaments allowing thin structure reconstruction [20], symmetric reconstruction [21] and even adaptive mesh refinement [22].Attempts have also been made to couple the level set method with the MOF method [23,24].Due to its expensive computational cost, further improvements have been made only for Cartesian cells using pre-computed values to find the reconstructed centroid in an efficient manner [25] and even a machine learning approach [26] improving drastically the runtime.When large deformations of the interface occur, standard MOF techniques are not precise enough to maintain a smooth interface and breakup occurs similar to VOF methods.To overcome the problem of unphysical breakup, we propose a novel symmetric multi-material approach to maintain the morphology of the interface for under-resolved filamentary structures.By combining the advantages of each approach, we construct a more precise interface at maximum deformation while maintaining an acceptable runtime.The optimisation algorithm uses a bisection method that does not require any parameter tuning.Capturing and reconstructing the exact topology yields time-consuming computation which has been simplified in our model.The novelty of the proposed method and the difference between the filament and standard MOF methods are highlighted in Fig. 1.A new test is presented to highlight large deformation of thin interfacial structures.
The paper is structured as follows.In Section 2, the standard MOF and its advection approach are described detailing the choice of using a bisection method in the optimisation part.Then, the filamentary method is presented in Section 3. The conglomeration algorithm and adjacency test are described, as well as the choice of capping the number of materials at three.Finally, Section 4 presents and analyses the results for many advection benchmark problems to demonstrate the accuracy and advantages of the proposed method.Some concluding remarks are given in Section 5.

Problem definition
Let us define the problem imposed by the MOF method in order to reconstruct an interface.Consider a convex polygon ω that is defined by n vertices, x 1 ,..., x n .The area of ω, denoted |ω|, and the centroid (centre of mass), denoted x c (ω) can be computed as follows (1) Note that x n+1 = x 1 .Let depict an arbitrary convex cell, hence not restricted to a Cartesian cell, filled with two different materials.Only considering the first material μ 1 , its area relative to the area of the cell is denoted by F ref (μ 1 ) which corresponds to the volume fraction.Similarly, x ref (μ 1 ) is defined to be the reference centroid of μ 1 within the cell.
The MOF reconstruction problem is formulated as an optimisation problem in which the distance between the reference centroid x ref (μ 1 ) and the reconstructed centroid x act (μ 1 ) is minimised while keeping the volume fraction of the reconstructed polygon F act (μ 1 ) equal to the volume fraction of μ 1 .One can summarise the optimisation problem as follows: If μ 1 already occupies a polygon with a piecewise linear interface, the MOF method aims to reconstruct the exact interface.As shown in Fig. 2, the reference interface may be curved, hence the minimised centroid distance will aim to give the best reconstruction.

Reconstruction
The reconstructed normal to the interface within a polygon can be evaluated analytically [17] but only for rectangular cells.However, for cells of any other geometrical shape, a minimisation algorithm is needed to evaluate the normal to the interface.The unit normal is defined to be n = [cos(φ), sin(φ)] where φ corresponds to the angle the interface makes with the horizontal.To cover all possible normals, φ ∈ [0, 2π ].The minimisation function, also known as the objective function, is recalled In general, f (φ) may have multiple local minima.The first derivative of the objective function for a convex cell, initially given in [16], is defined by where x act (φ) is given by and is evaluated using the length of the reconstructed interface segment (φ).

Bisection method
In this section, a new algorithm is presented to evaluate the normal to the interface.The algorithm used in this paper to find the global minimum is a bisection method.Using four quadrants, explicitly [0, π/2], [π /2, π], [π , 3π /4] and [3π /4, 2π ], the zeros of the first derivative of the objective function can be determined.The bisection method uses only a maximum of 10 iterations per quadrant to find the local minimum with a tolerance of 10 −10 .When the value of the first derivative falls below the specified tolerance at the boundaries of a quadrant, the bisection method is terminated for that quadrant.Once the minimum for each quadrant is found, evaluating the objective function for all valid values will give the global minimum.The global minimum of f (φ) will result in the best approximation for the optimisation problem defined above.Fig. 3 shows the set of solutions as well as the objective function within the four quadrants.Knowledge of the normal enables one to flood the cell [17] to reconstruct the interface with the minimum distance between the reference and reconstructed centroid, which is defined as the least centroid error.This method has the advantage of not requiring any initial condition and fine parameter tuning to converge to the solution and is guaranteed to find the global minimum.However it may require a larger number of iterations to converge.

Advection
Dynamic tests involve advecting materials across multiple time iterations.Information from the previous time step is needed in order to reconstruct the material interface at the next time step.The most natural way to perform this reconstruction is to use a Lagrangian pre-image to capture the volume fraction and centroid of a material.
All vertices of a cell are advected backwards in time using a 2 nd -order Runge-Kutta scheme to form the backtrace cell as seen in Fig. 4a.The backtrace cell may intersect several cells at the previous time level.The Sutherland-Hodgman polygon clipping algorithm is used in order to intersect each of these cells to gain information about volume fraction and centroid.The advantage of using the Lagrangian approach is that there is no limitation on the CFL number used in the model.Moreover, the Lagrangian advection procedure is said to be unsplit, which means it only requires one advection and reconstruction per cell [27].

Advection of volume fraction
To compute the volume fraction at the next time step, the sum of intersecting areas forms the new volume fraction of the cell as highlighted in purple in Fig. 4b.However, in some cases, its value may depend on the backtrace cell area relative to the cell area.If the backtrace cell area is larger than the cell area, there is potential for the volume fraction to exceed unity.On the contrary, if the backtrace cell area is smaller than the cell area, there is potential for the volume fraction to be  smaller than unity while being entirely filled with one material.These cases may occur when the backtrace cell intersects with only one material, making the new theoretical volume fraction equal to unity but the actual volume fraction is either greater than or less than unity.If this is the case, a post advection remapping procedure is introduced in order to ensure that the total material mass is consistent throughout the advection time.The difference between the actual volume fraction and unity are gathered, then redistributed equally across all cells that can accept a gain or loss of mass/volume fraction.This is defined as a global redistribution [28].The modified mass in each cell is negligible so that the shape of the interface is not changed significantly, which has been demonstrated later in the validation.To be less expensive, this procedure is only performed once per time step, which means there is a risk of not being able to redistribute the total mass.

Advection of centroid
To compute the centroid at the next time step, the centroid of the intersection of the backtrace cell with a cell is computed, then advected using the same scheme as for the backtrace cell advection as shown in Fig. 4c.All cell intersection centroids are advected forward in time individually.The new reference centroid is obtained by weighting the cell intersection centroids with their volume fraction.Since all centroids are framed within the backtrace cell at the previous time step, the new reference centroid is guaranteed to be within the cell after forward advection.

Filament capturing
Filaments are defined as thin strands of material surrounded by another material within a cell.These are structures thinner than a cell size.A standard MOF reconstruction would create a linear interface splitting the cell in two, hence not reconstructing the topology correctly as shown in Fig. 5.When considering a filament, two linear interfaces emerge, one on each side of the structure, meaning that two reconstructions are needed to capture the topology perfectly.In filament reconstruction, the conglomeration algorithm is capable of detecting polygons of the same material that are not adjacent by using the numerical adjacency condition.A fictitious material is introduced to reconstruct one of the polygons surrounding the filament.Once reconstructed, the fictitious material is reassigned to its original material.A symmetric multi-material reconstruction is presented to generate a better topology.

Conglomeration
Filament reconstruction is performed when some adjacent polygons forming one material, called a conglomerate, are not adjacent to other conglomerates of the same material.The conglomeration algorithm allows us to identify whether a cell needs a multi-material reconstruction or a standard reconstruction.It is possible to identify all polygons of one material intersecting with the backtrace cell as shown in Fig. 6.Once all of these polygons are gathered, the conglomeration algorithm tests if each of these polygons are adjacent to each other.Conglomerates are considered even when they do not split a cell, i.e. being only adjacent to one cell edge.The green conglomerate in Fig. 8b is one of these.Flotsam are not discussed in this paper, in general on a coarse mesh they do not tend to exist.If more than one conglomerate is found, then one of these conglomerates is considered to be the fictitious material.The conglomeration algorithm is a tree-based structure testing adjacency of a list of polygons until the lowest level does not find any adjacent polygons.Algorithm 1 details the procedure to identify conglomerates.The reference volume fraction and reference centroid can easily be computed.

Adjacency
The adjacency test is performed on all sides (segments) of a polygon with respect to another polygon.Some tolerance is accepted as sides may not be perfectly adjacent but can still be considered adjacent.For the purpose of numerical round-off errors, each segment is described by a vector and if the magnitude of the cross-product of two vectors meets the lower bound of a tolerance, here x y with = 10 −3 , then segments are considered parallel.Segments may be considered parallel, yet they also need to be adjacent.Hence, the endpoint of a segment is projected onto the line defined by the other segment.If the distance between the endpoint of the segment and its projection falls below the specified tolerance, here x , the projection of an endpoint also needs to fall between the bounds of the other segment.Only if all conditions are satisfied are both polygons considered to be adjacent and hence form a conglomerate.Algorithm 2 summarises the conditional procedure to test if two polygons are adjacent with three nested conditions.Fig. 7 shows two polygons within a cell.Segments are highlighted in order to indicate the process of evaluating parallel and adjacent segments from two distinct polygons.Condition 1 is represented with gold segments.Condition 1 and 2 are represented with blue segments.All three conditions are represented with red segments.

Limitation to three materials
It may happen that more than three conglomerates form within the backtrace cell.In that case, a multi-material reconstruction can be considered.However, it can lead to expensive reconstruction when testing all the combinations for several cells per iteration.For this purpose, the number of conglomerates is capped at three in our model.Conglomerates are sorted by their volume fraction.
If two conglomerates exist for each material, the following condition is tested.If the second conglomerate of one of the materials has a volume fraction smaller than 10 −3 , then its volume fraction is added to the main (largest in volume fraction) conglomerate.If there still exists two conglomerates for each material, no conglomerates are considered and a standard reconstruction with the total volume fraction per material is performed.Fig. 8a highlights this scenario.Indeed, coloured conglomerates belong to Material 2, explicitly 2 and 2 .Material 1 also has two conglomerates in white, explicitly 1 and 1 .None of them are smaller than 10 −3 in volume fraction.
In other cases, conglomerates with the smallest volume are "reattached" to the largest conglomerate of the same material in the cell, usually where one material has one conglomerate and the other has more than two conglomerates.Then, these smaller conglomerates have their volume fraction added to the largest conglomerates.Fig. 8b highlights this scenario.Three conglomerates (coloured) belong to Material 2, here explicitly 2, 2 and 2 .Conglomerate 2 will be reattached to Conglomerate 2, while Conglomerate 2 will be considered to be the fictitious material for reconstruction."Reattaching" to the nearest conglomerates based on the distance between their respective centroids may also be considered but does not affect the topology greatly as the volume fraction for these conglomerates is often very small.

Symmetric reconstruction of filaments
The reason to cap the number of materials at three is based on computational cost.Reconstructing more than three materials at once has a significantly higher cost than only three materials.In addition, using a symmetric reconstruction of filaments may provide a better topology in material reconstruction.
A standard reconstruction aims to reconstruct an interface based only on minimising the centroid error of one material regardless of the other material in cell reconstruction.In some cases, this can lead to a large error in the remaining material centroid.The symmetric reconstruction approach aims to minimise both centroids at the same time.The objective function f sym (n), combining both centroid defects, is given by where x ref (μ rem ) denotes the reference centroid of the remaining material in cell, and x act (μ rem ) is its reconstruction centroid.
When it comes to filament reconstruction or three material reconstruction, the standard approach is to test all ordering combinations and evaluate the topology that reduces the total centroid defect.This procedure is called a sequential reconstruction.The total centroid defect E can be expressed as the sum of all material μ i centroid errors Consider three materials A, B and C, then six different configurations are possible.Explicitly, and in order of reconstruction, these are (ABC), (ACB), (BAC), (BCA), (CAB) and (CBA).A symmetric reconstruction reduces the number of combinations to only three, thereby reducing the computational effort.Considering the same materials, (ABC) and (ACB) would be redundant as the first reconstruction minimises A and the grouping of B and C.Then, (BC) or (CB) will result in the same reconstruction as only symmetric reconstruction is considered.As seen in Fig. 9, a symmetric reconstruction provides a better topology.

Results
In this section, several benchmark problems are considered with the aim of testing the performance of the new filament MOF method.Several problems are of considerable interest since the associated velocity field yields high deformation in the material.Maintaining the correct topology at maximum deformation is attractive and important for most engineering problems.However, in order to assess the predictive capability of interface capturing methods, each of the flows is reversed over the same time period and compared to its original configuration.Whilst comparison with the initial condition is possible, the MOF enables one to evaluate the difference between the final reconstruction and the original/reference configuration rather than the initial reconstruction.From a computational cost perspective, our model uses an analytical reconstruction where possible.Indeed, when only the reconstruction of a piecewise linear interface between two materials in a Cartesian cell is required, this approach is significantly more efficient.In order to reduce the total error in reconstruction, the interface is reconstructed based on the material with the smallest volume fraction in a cell as suggested by Makundan et al. [24].For cases involving more than two materials reconstruction, the symmetric multi-material approach is chosen.

Error evaluation
Evaluating errors in interface reconstruction is a powerful tool to compare different interface tracking/capturing methods.The numerical errors in terms of volume fraction can be evaluated with the L 1 error norm E L 1 as its relative error norm E r and its maximum error norm A more representative error measure is the symmetric difference error which provides a better estimate of the interface reconstruction error.The symmetric difference error E sym is given by (11) When comparing the reference interface with its reconstruction in individual cells, the symmetric difference error can be interpreted as the area between the two interfaces.Fig. 10 shows three different scenarios of intersecting interfaces and highlights the area corresponding to the symmetric difference error E sym .Some simple calculations are necessary to evaluate the area of a segment.
As well as evaluating the error in reconstruction, ensuring mass conservation is also crucial during these advection tests.In 2D, mass conservation corresponds to area preservation and mass loss is given by the expression

Benchmark: Zalesak slotted disc
In this benchmark test case, a slotted disc is advected in a rigid body rotation motion.A circle of radius r = 0.15, with a slotted rectangle of width w = 0.05 and a maximum height of h = 0.85, is centred at (0.5, 0.75) in a unit domain.The corresponding velocity field is given by u(x, y) = 0.5 − y x − 0.5 (13) This case does not exhibit any filament formation, however it shows that the conglomeration algorithm works for velocity fields that are rotating rather than deforming the interface.Five different uniform grids have been used explicity 32 × 32, 64 × 64, 128 × 128, 256 × 256 and 512 × 512.On the coarsest mesh 32 × 32, the number of iterations is set to be n it = 1256 and t = 2π /n it .The number of iterations is increased proportionally with increasing mesh refinement.Hence, t is decreased correspondingly.
In order to study the error convergence, the error calculation for this test case is based on the L 1 error.Table 1 summarises the L 1 error for different mesh sizes and highlights the convergence of the numerical approximation.In addition, Fig. 11 highlights the final solution after a full body rotation for the first three grids.The shape of the original interface is captured well.However, the sharp edges around the slotted rectangle have been smoothed out during the rotation.Indeed, the MOF method is not able to capture these edges regardless of the degree of mesh resolution.The maximum error L ∞ is a more relevant measure of the error for this problem in order to understand the order of convergence around sharp edges.In this case, second order convergence may be attained in some instances but it may depend on the alignment of the sharp edges of the slotted disc with the grid.
The behaviour of the L 1 error over one rotation is shown in Fig. 12 for three different meshes mentioned in Fig. 11.The plot highlights that despite the interface only rotating, the error increases during the rotation progresses as the interface reconstruction error accumulates at each time step.

Benchmark: reversible vortex T = 8
The reversible vortex is a benchmark test case for deforming advection cases.A circle of radius r = 0.15 centred at [0.5, 0.75] in a unit domain is deformed in a divergence-free velocity field given by u(x, y, t) = − sin 2 (π x) sin(2π y) sin 2 (π y) sin(2π x) cos(πt/T ) (14) where T represents the full period and T /2 the time at maximum deformation.Here T = 8.The Courant-Friedrichs-Lewy (CFL) number is 1, hence the number of iterations n it = 256 and t = x when a 32 × 32 uniform Cartesian mesh is considered.The number of iterations increases proportionally with the mesh.
The circle deforms in a filamentary structure at maximum deformation t = T /2.For this test case, filament detection is enabled.Several grids from 32 × 32 to 1024 × 1024 have been used to perform this dynamic test case.

Table 2
Reversible vortex test case data using T = 8 compared with the standard MOF (STD MOF) and with results generated using other MOF methods in the literature: a standard MOF with adaptive mesh refinement (AMR) method [22], a filament AMR method [20], a coupled level-set MOF (CLSMOF) [24].The symmetric difference error, E sym , is shown in Table 2 for the initial reconstruction and at the final stage.The performance of our method is compared with the results obtained using other MOF methods as well as with the standard MOF.Runtime, rounded to the next integer value, is also compared because the MOF method can be computationally expensive.Currently, the code has not been parallelised and so the computations are performed on a single core.The order of convergence of this method is also highlighted as well as the mass difference.Fig. 13 shows the maximum deformation before reversal and the final reconstruction for different mesh sizes: 32 × 32, 64 ×64 and 128 ×128, respectively.Using a filament approach, the vortex does not exhibit any spurious separated structures, even on a coarse mesh.In this test case, the trailing tail shows a thicker structure as the coarse cell cannot reconstruct the filament tail accurately.As the mesh is refined, the tail becomes well-defined but thicker than the filament width.The MOF method naturally creates these structures as it exhibits some cross-stream diffusion, leading to a shorter tail than expected.The symmetric difference error converges slightly faster than other MOF methods with a smaller error on the finest mesh.Runtime is also considerably faster by a factor of between two to five.The error shows high order of convergence, almost matching the reference order two.The symmetric difference error for standard MOF exhibits a slower order of convergence on coarser grids.However, the symmetric difference error is almost indistinguishable on the finest meshes for the two approaches.

Influence of the mass redistribution
The remapping procedure does not affect the topology greatly as the volume fraction that needs to be redistributed during the procedure is very small.Indeed, during the deformation of the vortex the volume fraction redistributed oscillates between 10 −4 and 10 −10 .During the early stages of the deformation, most mass has to be redistributed as there are many cells in the inner part of the circle that are over/under-filled and very few cells are mixed cells, i.e. cells containing an interface.On the contrary, at maximum deformation, very few cells are over/under-filled cells, most of them contains an interface, or two in the case of filaments.Fig. 14 summarises this.The difference between the interface shapes obtained with and without the post advection remapping procedure is highlighted in Fig. 14a, whereas Fig. 14b shows the variation of volume fraction redistributed per iteration.Note that the mass is redistributed equally between mixed cells.Table 3 shows the mass difference for the reversible vortex case (T = 8) for a case where the mass was redistributed and when it was not.
Two orders of magnitude of difference can be observed, which shows the advantage of using the proposed method for mass conservation.

Influence of the CFL number on the interface
One expects the CFL number to influence the interface reconstruction.However, the Lagrangian advection procedure is not greatly affected by the CFL number.Therefore, most cases are performed with the maximum available CFL number which equals unity.In theory, a CFL number greater than unity can be used for such advection benchmarks.However, the stencil used in the dynamic test procedure encompasses only a 3 × 3 stencil and therefore limits larger CFL numbers.Fig. 15 shows a zoom on the final reconstruction for different CFL numbers 0.2, 0.4, 0.5, 0.8 and 1.0 and the error convergence.A lower CFL number will induce a larger number of iterations, therefore increasing the chances of error in reconstruction.However, the difference in error is relatively small in magnitude.The difference on the interface only occurs near the top of the circle which is near the tip of the filament at maximum deformation.

Influence of the filament capable method
Filament capable MOF is able to reconstruct a moving interface with a greater accuracy and better topology.Indeed, under strong deformations, materials tend to break up when they are not predicted to do so.At the instant of maximum deformation, a continuous interface is more likely and will result in better modelling of multiphase flows.Fig. 16 highlights both visual reconstruction and convergence of the standard and filament solution.Fig. 16a shows that several break ups of the dynamic interface occur when a standard MOF reconstruction is implemented.The final reconstruction does not match the reference circle.Fig. 16b compares the order of convergence between a standard and the proposed filament approach, together with other MOF methods.Note that for the finer grids, the error tends to the same values as the thickness of the structure is larger than a cell size, hence the filamentary approach is not used as frequently during the dynamic test.

Test case: reversible vortex T = 12
This is the same benchmark test case as considered in Section 4.3 except that the full period is increased to T = 12.
A larger period increases the deformation and thinner filaments are exhibited.Table 4 summarises the symmetric difference error for six different mesh sizes from 32 × 32 to 1024 × 1024.As the material is more deformed than in the previous benchmark with T = 8, the expected symmetric difference error is larger.As the mesh is refined, there are no longer significant benefits associated with using the filament method as the thickness of the deformed filament is greater than a cell width.Consequently, the order of convergence decreases from quadratic to linear until the order of convergence of the filament MOF follows that for standard MOF.The mass difference is very comparable.However, runtime is increased significantly.Indeed, the number of cells containing a filament structure compared to a standard interface is very large.Fig. 17 highlights the morphology of the very thin interface.Because filament reconstruction is computationally more expensive, the runtime is increased by a factor of three.

Benchmark: droplet flow
The droplet flow test case has a nonlinear divergence free velocity field.The deformation of material tears an initial circle of radius r = 0.125 centred in a unit domain into a V-shape.The velocity field is given by The velocity field is a function of time as the amplitude, f (t), varies in time according to At T max = 0.8, time at maximum deformation, the flow is reversed.The flow is reversed smoothly during a transition period of t = 0.1.
This test case provides a good insight into the filamentary formation of materials as the filament tip is leading as opposed to trailing in the previous benchmark.For the base grid, 32 × 32, the number of iterations is set to n it = 160 for the entire simulation and t = 0.01.The number of iterations is increased proportionally with the mesh and therefore t is decreased proportionally with the mesh.
The dynamic test is performed for different grids from 32 × 32 to 256 × 256 using a filamentary method.As the mesh is refined, this approach becomes less relevant.The symmetric difference error is compared with [20] despite an AMR capability being used in that paper.In addition, details of mass conversation and runtime are given in Table 5. Fig. 18 shows the maximum deformation and final reconstruction for 32 × 32, 64 × 64, 128 × 128 grids, respectively.It can be seen that coarser meshes lead to larger error in reconstruction.In addition, the method exhibits some diffusion in the sense of "floating" elements.These "floating" elements could be attenuated with a higher tolerance in available cell volume fraction.Lower volume fraction tends to create long and thin polygons, hence a larger error in reconstruction.The lower bound of volume fraction available in a cell is set to 10 −5 in our model, compared to 10 −8 in most comparative studies.We note that both maximum deformation and final reconstruction show a symmetric left-right deformation.As the tip of the filament gets thinner, even the filamentary approach cannot reconstruct the structure accurately.This leads to a shrinked filament structure.When the grid is refined, the tip of the filaments are well-defined and the final solution shows acceptable errors.The mass difference is acceptable, bearing in mind the choice of only one round of redistribution.In terms of runtime, the performance is compared with AMR [20] which uses fewer cells than in this paper.The order of convergence of the solution shows a remarkable performance compared to methods described in other papers.

Benchmark: rotating filament
The rotating filament benchmark is a test case where a thin rectangle is advected anti-clockwise in a rigid body rotation motion.The velocity field is the same as in Sec.4.2 and is given in Eq. ( 13).The rectangle is centred at (0.505, 0.75) in a unit domain.Its initial width is w = 0.006 and height is h = 0.3.For a coarse mesh, here a 100 × 100 grid, the initial condition may already contain a filament structure.The corresponding number of iterations is set to n it = 300 and t = 2π /n it .
This benchmark can only be tested with a filament enabled approach on such coarse meshes.Indeed, even with a 200 × 200 grid, the filament body is subject to under-resolved filamentary structures.Fig. 19 shows the rotating filament at different stages of the full body rotation.The filament body is well reconstructed.However, both ends of the filament show cross-stream diffusion because the MOF method cannot reconstruct sharp edges accurately.In addition, the filament height is shortened due to the reconstruction error.The filament shortening matches with the height shown in [20].The zoom on the top left of the figure highlights both shortening of the filament and cross-stream diffusion compared to the reference rectangle outlined in black.

Benchmark: S-shape
The S-Shape test case comprises a circle of radius r = 0.25 initialised in the centre of a unit domain.The associated velocity field is nonlinear and divergence free and given by u The advection process creates a highly deformed and thin structure which means the filamentary capability is also enabled here.The amplitude f (t) is given in Eq. ( 16).However, in this benchmark problem the maximum deformation occurs at T max = 4 and the smooth transition period is t = 2.
This case shows strong deformation and thin structures, mainly in the centre of the domain.A coarse mesh would struggle to reconstruct these structures.Indeed, for the 32 × 32 grid in Fig. 20, the central part may have three interfaces within a cell.Therefore, capping to three materials is a limiting factor, creating larger errors in reconstruction.Because of large reconstruction errors, several structures may merge and lead to different end results.The 64 × 64 grid is fine enough to have a maximum of two interfaces in a cell.The deformed interface shows an accurate representation at maximum deformation.

Conclusions
In this paper, a new MOF method with a symmetric multi-material approach has been presented where thin structures are resolved using a filament approach for a fixed coarse mesh.A novel robust approach to solve the optimisation problem is proposed using a bisection method.No initial condition or parameters are necessary and the global minimum can be always found.A Lagrangian backtracking approach ensures that there is no limitation on the CFL number when advecting materials.Solving under-resolved filaments inherently involves a higher computational cost, which is reduced by choosing to cap the number of conglomerates at three and using a symmetric approach.As a result, almost quadratic order of convergence is achieved and the error converges as the grid is refined.However for complex and large material deformations, the limitation of this method is shown and the topology might not be well maintained at sharp edges.
This efficient approach is applied to several benchmark problems with different levels of deformation.Most of these benchmark problems are compared with different MOF approaches, filaments, AMR, CLSMOF and standard MOF.First, the Zalesak slotted disc does not exhibit any filament behaviour, yet our approach shows good qualitative results.Other benchmarks such as the reversible vortex and the droplet flow case are tested for large deformation highlighting the quality of reconstruction of filaments on coarse meshes.Finally, the rotating filament benchmark is presented, only applicable using filament reconstruction for such coarse grids.The limitation of our method is shown in the S-shape deformation benchmark.For most benchmarks, the error and runtime are comparable or better than other MOF methods.Furthermore, the accuracy in interface reconstruction is improved for large deformation.In addition, runtime has been decreased compared to most MOF methods.
The MOF method, like most interface capturing methods, diffuses when advecting sharp edges.In addition, the tip of filaments is not well-resolved regardless of the mesh resolution.In future work we would like to include the possibility to reconstruct four conglomerates within a cell while maintaining a high level of accuracy in reconstruction with an acceptable runtime.For a fixed coarse mesh, this may lead to an increased precision in thin layered filaments while reducing the natural diffusion of material.This approach could involve an optimised selection of which material to reconstruct.We would also like to use a selective adaptive mesh refinement method for complex and large deformation.Coupling our accurate MOF method with a fluid flow solver is our next priority with interest in both incompressible [29] and compressible [30] multiphase flows.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Flowchart highlighting the key steps for a filament MOF method compared to a standard MOF method.

Fig. 2 .
Fig. 2. Reference vs. reconstructed interface with their respective centroids x ref and x act .n denotes the reconstructed normal to the interface.denotes the length of the interface segment.

Fig. 4 .
Fig. 4. Dynamic test: advection of backtrace cell backwards, intersection of volumes, advection of centroids individually.(For interpretation of the colours in the figure(s), the reader is referred to the web version of this article.)

Fig. 6 .Algorithm 2
Fig. 6.Conglomeration of polygons within the backtrace cell (dashed black outline) leading to the creation of a fictious material for a 3-material reconstruction.(a) Material μ 1 has 1 conglomerate (outline in blue); (b) Material μ 2 has 2 conglomerates (outline in green).

Fig. 7 .
Fig. 7. Schematic diagrams to test adjacent segments with another polygon: (a) shows two configurations where segments are parallel.Projection of the endpoints does not fall within tolerance (highlighted in gold).One of the projections of the endpoints does not fall within the other segment (highlighted in blue); (b) shows two segments that are parallel and adjacent, the projection of the endpoints falls within tolerance and within the other segment.

Fig. 8 .
Fig. 8. Schematic showing two complex examples of sorting multiple conglomerates within the same cell.We assume all coloured polygons belong to Material 2.

Fig. 9 .
Fig. 9. Comparison of (a) sequential and (b) symmetric reconstruction when using three materials, where (+) denotes the reference centroids and (o) denotes the reconstructed centroids.

Fig. 10 .
Fig. 10.Symmetric difference error E sym in a single cell.The area shaded in blue highlights the error corresponding to E sym irrespective of the number of times, n, the reconstructed interface intersects the reference interface: (a) n = 2, (b) n = 1, (c) n = 0.

512 4 .Fig. 11 .
Fig. 11.Solution of rigid body rotation for the Zalesak slotted disc.Green depicts a quarter of rotation.Blue half rotation.Purple three quarter of rotation.Red depicts a full rotation and final solution.The black outline depicts the reference interface.

Fig. 12 .
Fig. 12. Behaviour of the L 1 error during the rigid body rotation of the Zalesak slotted disc for different mesh sizes.denotes the angle of full body rotation.

Fig. 13 .
Fig. 13.Reversible vortex test case using T = 8 for 32 × 32, 64 × 64 and 128 × 128 grids.Top row of figures shows the maximum deformation.Bottom row of figures shows the final interface.

Fig. 14 .
Fig. 14.Comparison showing the effect on the interface shape of the post advection remapping procedure for mass conservation and the actual mass redistributed per iteration for different grids.

Fig. 15 .
Fig. 15.(a) Influence of the CFL number on the final reconstruction of part of the interface; (b) symmetric difference error, E sym , as a function of the CFL number for two grids: 32 × 32 and 64 × 64.

Fig. 16 .
Fig. 16.(a) Influence of the filament capable method on the reconstruction; (b) Symmetric difference error E sym compared with other MOF methods.Convergence rate is compared with a linear and quadratic reference.

Fig. 18 .
Fig. 18.Intermediate and final reconstruction for the droplet flow test case for different mesh sizes.Red depicts the maximum deformation before reversal.Green depicts the final reconstruction.The black outline is the reference circle.

Fig. 19 .
Fig. 19.Solution of rigid body rotation for the rotating filament.Green depicts a quarter of rotation.Blue half rotation.Purple three quarter of rotation.Red depicts a full rotation and final solution.The black outline depicts the reference interface.

Table 1
Dependence of the L 1 error, E L1 , relative error, E r , and maximum error L ∞ on mesh size for the Zalesak slotted disc problem.

Table 3
Mass difference for the reversible vortex with and without post advection remapping procedure.

Table 5
Symmetric difference error, order of convergence, mass difference and runtime for the droplet flow test case at final reconstruction compared to reference papers.