Visual analysis of interface deformation in multiphase flow

In multiphase flows, the evolution of fluid-fluid interfaces is of interest in many applications. In addition to fluid dynamic forces governing the flow in the entire volume, surface tension determines droplet interfaces. Here, the analysis of interface kinematics can help in the investigation of interface deformation and the identification of potential breakups. To this end, we developed a visualization technique using metric and shape tensors to analyze interface stretching and bending. For interface stretching, we employ the eigenpairs of the metric tensor defined for the deformation rate of the fluid surface. For interface bending, we present a technique that locally captures the interface curvature change in terms of a shape tensor, extracting its principal directions and curvatures. We then visualize interface deformation by combining both representations into a novel glyph design. We apply our method to study multiphase flow simulations with particular emphasis on interface effects. These include the interplay between fluid dynamics and surface tension forces leading to breakup processes following droplet collisions, as well as droplet-droplet interactions of different fluids where Marangoni convection along the surface is explicitly taken into account.


Introduction
A multitude of phenomena in science and engineering involve multiphase flows where a liquid phase is dispersed into another surrounding phase.The applications of liquid-gas configurations range from spray cooling systems in the food or microelectronic industry to droplet evaporation and splashing in combustion engines.Examples of liquid-liquid flow employment include waste and pollutant treatment in process engineering, transport of agents in medicine, as well as generation of powders and encapsulated material in the food industry.Due to their ubiquity, dispersed phases and their physics attract much research that focuses, among others, on the analysis of interface deformation, which can help to understand the relations between micro-scale and macro-scale phenomena.For instance, the ratio of the droplet's surface area to its volume affects the effectiveness of combustion engines and the amount of pollutant emissions, as a larger surface area results in faster fuel evaporation.Therefore, characteristics of droplet breakup and splashing on walls concerning atomization are essential research subjects.
The deformation of the interface between the liquid and the gaseous phase is caused by the dynamic motion of the involved fluids and by forces and stresses acting exclusively on the free surface of the liquid phase.The latter is caused by surface tension forces or by their components normal and tangential to the interface, respectively.The normal component of the surface tension force is proportional to the curvature of the interface and acts to minimize the surface area.Thus, it has a stabilizing effect on the interface as it acts against small surface undulations.If a liquid ligament is strongly deformed, this can also lead-in interplay with fluid dynamics forces-to a disintegration process, as the minimal surface area of a liquid ligament corresponds to the shape of a perfect sphere.In general, the surface tension coefficient depends on the temperature as well as on the chemical composition.A spatial variation of temperature or species concentration along the interface results in tangential stress.This stress induces a motion of the liquid at the interface, which is called Marangoni convection.Depending on the material, forces, and stress due to surface tension and fluid dynamics forces have different influences on the overall interface kinetics.
To understand phenomena arising at the interfaces, analysis of the aforementioned forces is particularly important, and visualization lends itself well to investigating their interdependence.However, visualization of multiphase flows is challenging due to complex dynamics, including collisions and the interplay of different densities and viscosities at the fluid interface.Moreover, the analysis of interface deformation requires an approach that can capture its important characteristics and effectively present them.Therefore, we have developed a visualization technique based on metric and shape tensor analysis that allows for efficient investigation.On the one hand, the metric tensor (Section 4.1) provides information on surface stretching and contraction, employing their respective directions and corresponding extensiveness, and can highlight possible breakup regions.This way, we can analyze the influence of surface tension force on droplet cohesion and breakup (Sections 5.1 and 5.2).The shape tensor (Section 4.3), on the other hand, indicates shape deformation in terms of curvature change and is used to provide insights into interface kinetics.Thus, it allows, e.g., domain scientists to observe the propagation of a capillary wave for the analysis of Marangoni convection (Section 5.3).For simultaneous visualization of interfaces stretching and bending-in short, interface deformation-we introduce a specialized glyph design (Section 4.4).The main contributions of our paper include the following: -computation of interface bending based on the shape tensor, -design and visualization of an interface deformation glyph, -analysis of the influence of surface tension on phase instabilities, and -application of our method for the analysis of the propagation of capillary waves caused by Marangoni convection.
Please note that in multiphase flow, interface stretching does not have the same physical meaning as in solid materials.Whenever the surface area increases due to deformation, inner fluid molecules move to the outer region, contributing to the additional interface area.Nevertheless, such a formulation allows for a practical description of interface deformation and is suitable for analysis on the macro-scale.

Related work
The topic of material interface reconstruction in multiphase flows is strongly related to visualization.While algorithms may be able to reconstruct arbitrarily many interfaces within a single cell (Bonnell et al. 2003), the volume fractions enclosed within the reconstructed interfaces are generally not preserved.The marching cubes algorithm (Lorensen and Cline 1987), widely used for visualization, suffers from the same issue.This problem was addressed by introducing more accurate reconstruction techniques (Anderson et al. 2010;Meredith and Childs 2010) and analyzing the reconstructed interfaces' stability by comparing them with time surfaces (Obermaier et al. 2012).Further, the PLIC reconstruction (cf.Section 3) was visually analyzed, providing several measures to assess the quality of the approximation (Karch et al. 2013).
Several solutions have been proposed to determine the curvature of the liquid surface.Some of them are able to use scalar fields directly (Cummins et al. 2005) (e.g., VOF field, cf.Section 3).For this, they apply a function to the scalar VOF field, such as a height function, reconstructed distance function, or convolved VOF function.For other methods, on the other hand, the interface needs to be reconstructed in order to combine the height function approach with paraboloid-fitting to accurately and robustly approximate the curvature (Popinet 2009).
Surfaces and volumes have been extensively used in flow visualization for the analysis of deformations, e.g., employing glyphs and stream volumes to visualize deformations induced by vector fields (Obermaier et al. 2009).Metric tensors for glyph-based visualization of deformation can be used to analyze stretching and shearing (Obermaier and Joy 2012).For visualizing integral surfaces in divergence-free vector fields, stretching information helps adaptively find stretch-minimizing stream surfaces (Barton ˇet al. 2015).Metric deformation analysis parametrizes and visually compares time surfaces from different simulation time steps (Brambilla et al. 2016).However, these works address single-phase flow where the surfaces generally follow the vector field.In our work, on the other hand, the fluid interfaces undergo deformation not only caused by advection but also by phase dynamics and surface tension and therefore require a different approach for constructive analysis.For general information on flow visualization, we refer the reader to surveys on this topic (Laramee et al. 2004;Lips ¸a et al. 2012;McLoughlin et al. 2010).
In our work, we utilize the notions of surface parametrization and related metric and shape tensors for the description of deformation (Floater and Hormann 2005) and visualization of deformation (Berres 2015).The concept of curvature difference was utilized for the computation of discrete shells used in computer graphics and animation for modeling deforming thin objects (Grinspun et al. 2003).Especially glyphs lend themselves well to the visualization of tensor fields (Bi et al. 2019;Hergl et al. 2021), in our case for symmetric second-order tensors.One prominent example is superquadrics (Kindlmann 2004;Schultz and Kindlmann 2010), which combine the advantages of ellipsoids and cuboids into a shape that indicates the certainty of orientation.Because we provide two tensors, we decided to create a combined glyph, similar to the idea of probes (de Leeuw and van Wijk 1993).
We exemplify the utility of our technique using multiphase flow simulation data.The investigation of droplets is an active area of research (Focke and Bothe 2012;Heinemann et al. 2021;Jadidi et al. 2019;Tryggvason et al. 2011;Yokoi 2008), as is the analysis of liquid jets (Ertl and Weigand 2015;Fuster et al 2009;Li et al. 2010).We refer the reader to Fuster et al. (2009) for a detailed introduction to multiphase flow simulations and Lefebvre (1989) for a thorough description of liquid atomization and sprays.
This work started as a master thesis (Straub 2016) and was briefly discussed in project reports (Straub et al. 2018(Straub et al. , 2019) ) of a collaborative interdisciplinary project on droplet interfaces (International research training group 2160), as well as in an overview article (Straub and Ertl 2020).However, these publications only gave a broad overview of our technique and showed parts of the results.Here, our approach is described in more detail, and the results are thoroughly discussed.Additionally, this paper introduces a new visualization using specialized glyphs for interface stretching and bending.Please note that part of this work has been used in a dissertation (Karch 2017).
3 Multiphase flow data Hirt and Nichols (1981) introduced the volume of fluid (VOF) method for two-phase flows, where an additional scalar field is introduced to track the fluid interface between the phases.The initial approach suffered from numerical diffusion, and much effort has since been put into reducing it.Among today's most commonly used approaches is the piecewise linear interface calculation (PLIC) (Youngs 1984), which approximates the interface by planar patches (Rider and Kothe 1998).
The datasets used in this work were simulated with the computational fluid dynamics (CFD) solver Free Surface 3D (FS3D) (Eisenschmidt et al. 2016) for incompressible multiphase flow.This solver performs direct numerical simulations (DNS) on very fine rectilinear grids.In particular, it can accurately simulate capillary waves (Steigerwald et al. 2021), which can be induced by Marangoni convection in our application case (cf.Section 5.3).For interface tracking, the solver utilizes the VOF method, which stores an additional scalar field f (x, t) containing the volume fraction of the liquid per cell.This means a cell contains only the liquid phase when f = 1, only the surrounding phase when f = 0, and the interface for 0 \ f \ 1. Both, this scalar field and the cell-based velocity field u(x, t), are stored on rectilinear grids.To maintain sharp interfaces, PLIC is used to approximate the interface using piecewise planar patches with normals parallel to the gradient of f .In our approach, we use PLIC reconstruction to visualize the liquid surface.Thereby, we stay as close to the simulation as possible and avoid topological inconsistency that would occur with other methods, such as marching cubes.

Visualization method
Various forces acting on the fluid phase cause changes to the interface shape.Surface tension induces spherical shapes, while fluid dynamic forces generally have the opposite effect often destabilizing the shape.
Phase deformation typically changes the interface area, which can be locally interpreted as interface stretching.Obermaier and Joy (2012) employed metric tensors to enhance surface flow visualization by providing stretching information.For the analysis of multiphase flows in our work, we also utilize the metric tensor.However, since the surface of interest is naturally defined by the fluid interface, the metric tensor has to be limited to the surface space.We describe the details of our approach in Section 4.1.
While the visualization of the metric tensor gives insights into the local variations of the surface area, the visual analysis of curvature change allows for a direct representation of shape deformation.It can be expressed as bending along two principal directions.Our approach extends the mean curvature change employed by Grinspun et al. (2003), who compute bending energy in discrete shells.Here, we additionally utilize the information on the directions of bending and their corresponding strength to provide more details on interface kinetics, as described in Sections 4.2 and 4.3.
For the visualization of stretching and curvature change, we introduce a specialized glyph.For both stretching and bending, this glyph conveys the respective principal directions and their corresponding magnitude.For this, the glyph consists of two separate parts, allowing for an intuitive interpretation of both tensors.Details are provided in Section 4.4.
Please note that in all our methods, we introduce for each cell a PLIC-local coordinate system P, its origin the respective interface barycenter x c .This system is spanned by the orthonormal basis e r ; e s ; e h h i , where e h is the PLIC plane normal, and e r , e s are chosen accordingly.

Metric tensor for stretching visualization
The interface deformation due to the simulation velocity u can intuitively be described by a metric tensor representing the first fundamental form I p from differential geometry in PLIC space P. For this, the deformation induced by the flow is approximated by the Jacobian J u [ R 393 of u.For two neighboring points P i and P j , the change in their displacement Dx = P j -P i during a single simulation time step Dt can then be computed as The matrix F = (I ?M) [ R 393 is the deformation gradient tensor.We can project the displacements Dx to the PLIC interface by right-multiplying F with N = (e r , e s ) [ R 392 .The resulting matrix J p = FN [ R 392 is the Jacobian of the displacement and, as noted by Floater and Hormann (2005), is used to construct the metric tensor as give the stretching factors for the corresponding eigenvectors.Hence, r 1,2 \ 1 indicates contraction and r 1,2 [ 1 indicates stretching in the respective direction.Finally, the stretching directions e 1,2 [ R 2 obtained from the eigenvectors of I p can be transformed to the original space by e 0 = Ne 1,2 [ R 3 , and-together with r 1,2 -can be used for visualization.Please note that since I p is a symmetric tensor, it does not capture deformation due to shearing.Also note that the component gradients in J u are computed using only points at the interface or within the fluid phase, which means that the Jacobian is estimated on scattered data.Therefore, we employ a regression-based method that accounts for unstructured point positions and use inverse distance weighting for improved accuracy (Correa et al. 2011).

Paraboloid-based surface approximation
For curvature computation, the approach described by Popinet ( 2009) uses the interface barycenters of the PLIC interfaces-or consistent interface positions calculated with the height function approach-from the central (x c ) and the neighboring cells (x i [ N x c ) to fit a paraboloid of the form pðr; sÞ ¼ ðr; s; hðr; s; fa i gÞÞ ð2Þ with the height function that depends on the coefficient set {a i } determined with linear least squares and is defined in the PLIC-local coordinate system P.It is worth noting that with the definition in Eq. 3, the vertex p(0, 0) of the paraboloid coincides with x c , irrespective of the coefficients {a i }.This is an approximation, as a general ansatz for h would contain additional terms a 3 r, a 4 s, and a 5 .Including a 5 would lead to different origins of the two paraboloids.Further, in the investigated datasets, the coefficients a 3 r and a 4 s had marginal effects on the computed quantities and only at interfaces whose least squares approximation was of low quality (e.g., for tiny droplets), so we neglect these terms in our further discussion, as well.
To analyze curvature change for a point x c , we fit two paraboloids p 1 and p 2 .First p 1 to the original point x c and its neighboring interface points x i [ N x c (Fig. 1, and green points in Fig. 2b), and then p 2 to the advected point set x i 0 (red points in Fig. 2c) as Here, we use the advected point set instead of the interface points of the subsequent time step because the point x c t does generally not correspond to a barycenter x c t?Dt of the interface at another time.The advection velocity is the fluid velocity in the simulation minus the velocity at x c and with possible rotations of the interface removed.Here, xi is the position relative to x c , defined as xi ¼ x i À x c .In Eq. 5, is the average angular velocity.With this procedure, we ensure that the two paraboloids p 1 (Fig. 2d) and p 2 (Fig. 2e) have the same orientation and share the same vertex x 0 c = x c .Hence, the difference of the paraboloids p d (Fig. 2f) is defined as p d = p 2 -p 1 = (r, s, h 2h 1 ).By computing the corresponding shape tensor, it can be used directly to extract information on the curvature change in terms of principal curvatures and directions, which correspond to bending strength and direction.The shape tensor is defined as , where the elements of the tensor I p are dot products of the partial derivatives of the surface function p(r, s) in Eq. 2, i.e., of the difference paraboloid p d , and the elements of II p are dot products of respective second-order partial derivatives and the interface normal e h , (II p ) ij = (q ij p) Á e h .For our purpose, we are only interested in S evaluated at x c .The first partial derivatives in r and s for the height function in Eq. 3 hence both vanish, so q i p Á q j p = d ij with Kronecker delta d ij , and I p is the identity matrix.The second derivatives read so, in summary, The Gaussian curvature K and mean curvature H are subsequently computed (Patrikalakis and Maekawa 2002) as From this, the principal curvatures, or eigenvalues of S, are obtained as The corresponding eigenvectors read k 1,2 = (a 2 , j 1,2 -2a 0 ) > .Please note that k 1 \ k 2 in all cases.The case a 2 = 0 has to be treated separately and gives j 1 = 2a 0 , j 2 = 2a 1 and k 1 = (1, 0) > , k 2 = (0, 1) > .

Interface deformation glyph
To visualize interface deformation, we combine two glyphs for interface stretching and bending, similar to the general idea by de Leeuw and van Wijk (1993).To create these glyphs, we exploit the tensors' properties, as well as the methods used to compute them.We discuss our design choices following the Design Guidelines (DG) stated by Borgo et al. (2013).
According to DG1, our visualization space is three-dimensional, although our glyphs are positioned along the surface of the fluid phase.Hence, the density of glyphs (DG2) and glyph placement (DG11) are constrained, as every cell containing the interface should contain a glyph.This also hinders the use of interactive methods, such as slicing and brushing (DG13).The fluid surface is always shown for context, leading to a hybrid visualization (DG3).As both glyphs encode orientation for stretching and bending, they are placed parallel to the surface and cannot benefit from billboard effects, as suggested in DG7.
For stretching, we can directly extract the two main directions (eigenvectors) together with their respective magnitudes (corresponding eigenvalue).Thus, we visualize stretching using superquadric glyphs.These glyphs are generally simple and symmetric (DG8), and mapping stretching to scaled superquadrics provides an intuitive (DG10) and perceptually uniform mapping (DG4).As our glyph lives in PLIC space, i.e., it is two-dimensional, we can reduce the explicit formulation for the glyph geometry by Kindlmann (2004) to two dimensions: where h [ [0, 2p] and a = (1c l ) c , with linear anisotropy c l = r 1 Àr 2 j j r 1 þr 2 , edge sharpness parameter c, and stretching factors r 1,2 with r 1 [ r 2 .However, we use this equation to form a disc glyph to allow the inner space to be occupied by our second glyph.As the glyph is embedded in three-dimensional space, we additionally add height to allow lighting to improve the perception of the glyph from different perspectives.This is in line with DG12 for enhancing depth perception and partly mitigates the issue of perspective projection when mapping to glyph size (DG14).The resulting glyph is shown in red in Fig. 3a.
Instead of the shape tensor, we employ the polynomial that is used to deduce said tensor to derive our bending glyph.This polynomial is the difference between the two polynomials fitted to the (advected) interface positions.Starting from a mesh representing a two-dimensional disk, we use the height function h p d (r, s) of p d to map the vertices (x, y) to (x, y, h p d (x, y)).This results in an intuitive (DG10), simple and symmetric glyph (DG8), whose height map h p d (r, s) provides a perceptually uniform representation (DG4).If there were no bending, we would therefore be left with a flat disk.Otherwise, a shape of would indicate increasing concavity, while would indicate increasing convexity in the respective principal direction.Again, lighting is paramount for depth perception of the three-dimensional glyph shape (DG12).Our bending glyph is depicted in green in Fig. 3a.
Because until now we map our two tensors exclusively to shape for attentive visual processing according to Borgo et al. (2013), we can further exploit other channels to guide the user's attention.According to studies on glyph design (Healey and Enns 2012;Ropinski et al. 2011), pre-attentive search-or a so-called pop-out effect-is best achieved with color.Therefore we map a value for importance to both color and opacity.This follows DG6 and DG5 by mapping our tensors not only to shape but redundantly also to color.To prevent colors from interfering with each other, we fix hue and value for both glyphs individually and only map the magnitude of stretching or bending to saturation in HSV color space, as well as to opacity.This provides necessary orthogonality and allows us to differentiate both glyphs intuitively (DG9).Note that our implementation is a ParaView plugin, and ParaView does not support HSV color maps.Hence, we approximate this behavior in RGB space.With this mapping, the user is able to quickly identify regions with either weak or strong deformation.As this mapping is done separately for both glyphs, it is also possible to differentiate regions where only one of the two deformation types dominates.For an example droplet, the visualization can be observed in Fig. 3b and for a zoom-in view in Fig. 3c.Because we have two magnitudes (eigenvalues) for both stretching (r 1,2 ) and bending (j 1,2 ), we have to combine those into a single value, respectively.For bending, we can trivially use the sum of the absolute changes in curvature ĵ = |j 1 | ?|j 2 |.However, for stretching, we have to multiply the maxima of the magnitudes and their inverse to get the combined magnitude for stretching rˆ= max {r 1 , r 1 -1 }Á max {r 2 , r 2 -1 }.This can be understood analogously to the absolute operator | Á | for summation, but now values centered around 1 are mapped from (0, 1] to [1, ?).
To allow the user to perceive stretching or contraction from glyph shape easily, we additionally provide an optional reference circle (Fig. 4a) in the visualization.With this, even weak stretching and contraction become distinguishable.An example application is shown in Fig. 8b.Further, both glyphs-for stretching and bending-can be extended to show their main directions using strips (Fig. 4b, c) along the glyph surface.This is helpful if the user wants to compare the directions of neighboring glyphs explicitly.Additionally, color mapping can be added for aid when glyphs feature only weak stretching or bending, and thus shape alone gives only a poor indication.

Implementation
We have implemented our visualization approach as a set of ParaView (Ayachit 2015) filters for stretching analysis with the metric tensor, curvature change with the shape tensor, and glyph visualization.In our case, using ParaView allows the domain experts to use our method in a familiar framework.While the implemented glyph visualization provides explicit geometry for all cells, our implementation allows for easy adaptation to the GPU.Instantiation of the glyphs would then be done for each cell by providing transformation matrices.Additionally, all methods are parallelized using OpenMP and MPI and can thus be easily used on supercomputers to analyze large datasets.Our implementation is integrated in the TPF framework: https://github.com/UniStuttgart-VISUS/tpf.

Results
In this chapter, we present the utility of our visualization approach on three multiphase datasets.The first two are collisions of water droplets in a gaseous surrounding.We use them to show that our method can be used to spot droplet breakup.In our third dataset, the onset of coalescence of a water droplet and an ethanol droplet has been simulated, taking concentration-induced Marangoni convection into account.Here, we highlight the differences between simulations where Marangoni convection was enabled and where it was neglected.Contrary to the first dataset, the goal here is to aid the domain experts in implementing the Marangoni term with our visualization method.

Colliding drops dataset
Our first dataset is a simulation of two water droplets in a gaseous surrounding, stored on a Cartesian grid (256 3 cells).An overview of this data can be seen in Fig. 5a-f.Initially, the two droplets move toward each other, having opposing velocities, and eventually clash in an off-center collision.The result is the formation of a disc-shaped structure, which successively disintegrates into many smaller droplets moving away from the center of impact.In Fig. 6, the visualization of the metric tensor for stretching and contraction is depicted.Note that we first focus on the effect of interface stretching and therefore only show the stretching part of the glyph.Color and opacity mapping reveals areas of strong stretching and contraction.Especially along the ligaments connecting the larger drops (zoom-in), strongly elongated glyphs can be observed.Perpendicular to the ligaments, i.e., in the vertical direction, the glyphs show contraction.This indicates that the ligaments are subject to thinning.On the other hand, glyphs are stretched along the tunnel-like structures, i.e., in the horizontal direction.Both observations suggest that the ligaments are collapsing, eventually resulting in disconnected droplets.This observation is confirmed by the time series shown below.Here, it is clearly visible that a breakup between the larger drops occurs.Furthermore, an important point is that, in this scenario, the glyphs help to clarify the actual role of the surface tension force in this observed development.As mentioned before, the surface tension force is proportional to the curvature of the interface.Thus, the forces are strong in the radial direction of the stretched ligaments and at the transition between the stretched ligaments and the larger drops.Especially the forces at the latter locations lead to a continuous flow of liquid into the larger bubbles, which results in a thinning of the stretched ligaments.As the glyphs show exactly this behavior, they help the domain experts by directly confirming that the observed deformation is fully surface tension driven.

Raindrops dataset
The second dataset again shows two water droplets in a gaseous environment in a Cartesian grid of 358 9 256 9 256 cells.This time, however, they move in the same direction but at different velocity magnitudes.
Once the larger droplet catches up with the smaller one, they merge and later disintegrate into a large amount of smaller droplets.The time series of the visualized raindrops is depicted as an overview in Fig. 5g-l.In Fig. 7, a simple structure from the Raindrops dataset is shown.In this case, glyphs are contracted predominately horizontally, which indicates cohesion between the left and the right part.This cohesion prevents the droplet from disintegrating.Vertically though, glyphs are stretched, which leads to elongation and, thus, a thickening of the tunnel-like structure.This is the opposite behavior compared to Fig. 6, described in Section 5.1.Once more, an overview of the subsequent time steps is given in the images below.Taking additionally into account the change in curvature, as shown in Fig. 3c, one can observe an increase in convexity horizontally and an increase in concavity vertically.Both the increase in convexity and concavity Fig. 6 Stretch tensor in the Colliding Drops dataset for a time step between Fig. 5d and e. Color and opacity guide the user to regions of strong stretching or contraction.Especially along the ligaments connecting the larger droplets (zoom-in), the tunnellike structure is being stretched (horizontal), while perpendicularly, they are being contracted (vertical).The time series below depicts the subsequent time steps fit the indicated process of volumetric growth of the center region.The strong curvature in this region again indicates that this development is fully surface tension driven, as confirmed by the glyphs.
A more complex example and droplet structure can be seen in Fig. 8. Here, three different areas are marked where different behavior can be observed.One region shows strong deformation, and the other two exhibit breakup.The region marked as À is being stretched vertically and contracted in the perpendicular, horizontal direction.This leads to a deformation of the main body of the structure, with a vertical elongation and a horizontal shortening.In the other two regions, `and ´, the contraction orthogonal and the stretching alongside the ligaments lead again to breakups.This can also be observed in the images below, depicting the following time steps.The zoom-in on region `in Fig. 8b additionally features reference circles to aid the user in better identifying stretching and contraction from glyph shape.

Marangoni convection datasets
Two simulations were performed to show how a variable surface tension changes the dynamics of a free surface and how the previously presented visualization approach can help to explain observed phenomena.In these simulations, a water droplet and an ethanol droplet of equal size are at rest and touching each other (Steigerwald 2023).The results of both simulations are depicted in Fig. 9, and comprise 512 9 256 9 256 cells in a Cartesian grid.Fig. 9a shows the initial state for both simulations, whereas Fig. 9b-e show the onset of coalescence.In the lower half of the images, concentration-induced Marangoni convection is neglected; in the upper half, it is taken into account.Without Marangoni convection, the coalescence is almost symmetrical and similar to the coalescence of identical liquids.Only a slight asymmetry is visible due to the different properties of water and ethanol.By taking into account the concentration-induced Marangoni convection, the dynamics of the free surface change significantly.
Water and ethanol have very different surface tension coefficients: Due to this, a strong surface tension gradient causes tangential stress acting towards the higher surface tension, i.e., towards the water droplet.As a consequence of the tangential stress, a thin ethanol layer spreads over the water droplet.This is accompanied by the generation of capillary waves propagating over the water droplet, a phenomenon observed in experiments by Thoroddsen et al. (2007).However, the surface's actual stretching and bending behavior during the onset of coalescence remains unclear.By applying the here presented visualization approach, the surface behavior can be analyzed in more detail.In the case of the simulation without Marangoni convection, the neck between the two droplets grows almost symmetrically.This is reflected in both stretching and bending: in the vertical direction stretching, and a decrease in convexity can be observed.In the horizontal direction, there is an increase in convexity in the center, and to the sides at the base of the connection, there is an increase in concavity.In the case of (I),d Marangoni convection, there is less growth but an additional region on the left droplet exhibiting strong deformation.The white arrow in d indicates the direction of this capillary wave, as also indicated for a zoom-in on (I) c in Fig. 11d The deformation of the surface is visualized in Fig. 10.In the bottom row (II) and in e, Marangoni convection is neglected, which leads to an almost symmetrical growth as described above.This is supported by the deformation glyphs, where the neck between the droplets exhibits stretching in the vertical direction (going around the neck), as well as a decrease in convexity.Together with an increase in convexity in the orthogonal direction near the center and an increase in concavity to the side, this clearly indicates the growth of the neck.The same behavior, to a far less extent, can be observed in later time steps for the Marangoni dataset in Fig. 10(I) and d.Additionally, a region on the left droplet exhibits strong deformation.This region of strong deformation is caused by a capillary wave that propagates from the neck between the droplets towards the left droplet.Fig. 11 shows this wave for different simulation time steps.One can separate this region into four areas, from À in front of the wave, over `a rising and ´a falling slope, to the wake of the wave.This is most prominent in the last time step in Fig. 11d.In the earlier time steps (a-b), the wake of the wave is too close to the neck region, where deformation is predominately caused by coalescence.Interestingly, the wave causes stretching and contraction, but clearly observable only for a narrow part of the wave.By adding strips for the major directions, we can see that these consistently align with the direction of the wave, even when deformation is relatively weak, e.g., stretching in `.By adding color to the strips for interface stretching, it is easier to observe that there is a slight contraction along the direction of the wave in this area.Additionally, in front of the wave, there is a part where stretching and contraction seem to diminish over time.This might be because the effective radius of the radially propagating wavefront around the droplet is increasing, thus dampening the wave.The interface bending paints a clearer picture.As can be seen in all time steps and annotated in Fig. 11d, in front of the wave, we exhibit an increase in concavity À, followed by an increase in convexity `, an increase in concavity ´, and finally again an increase in convexity ˆ.This behavior is expected of a propagating wave, as seen in Fig. 12.However, our method is especially helpful at the beginning of the wave formation when the wave itself is not yet visible but is captured by our interface bending visualization.

Discussion
This section deals with the experience and actual application of the proposed visualization technique by domain experts from the scientific multiphase flow community.The section can be seen as critical feedback, and it shall help readers better understand the viewpoint and motives of domain experts to use the proposed visualization technique.
First, it is essential to mention that actual investigated multiphase flows differ very often significantly from the simplified application scenarios presented in this work, which were kept simple for demonstration purposes.In many cases, the liquid phase is, for example, not modeled as a pure Newtonian fluid like water, for which most domain experts have developed an intuitive understanding of the flow dynamics.The liquid phase very often exhibits a complex rheological material behavior, which can even include a complicated time dependency, like in the case of thixotropic materials.Therefore, analyzing the deforming interface of such liquids is generally difficult.Especially for such complex liquids, the proposed visualization technique has proved to be an effective tool concerning the analysis of the often peculiar development of the interface.In these scenarios, the glyphs help to gain a deeper understanding of the liquid motion as the obtained interface deformation can be associated more easily with the spatially varying physical properties of the liquid phase.
At this point, it is also important to mention that for typically highly resolved production runs on HPC systems, the disk quota often limits the amount of data that can be written out during a simulation.A high temporal resolution of disintegration processes, like the ones presented in this work, tends to be the exception while also simulating at a high spatial resolution.
In our opinion, especially for these spatially highly resolved cases, for which the amount of output data is limited, and single events in the simulation are probably only shown in a single time step, the visualization technique displays its strength.The understanding of the superquadric glyphs for the directional increase of the interface, together with the brightness for the corresponding intensity, is very intuitive and straightforward for the dense information content of the glyphs.The glyphs that represent the curvature change are Fig. 12 Propagating capillary wave with sketched regions based on interface bending.The wave is moving to the left, where the current wave is drawn in black and the future wave in gray.À,´Increase in concavity, `,ˆincrease in convexity also very helpful but are not as intuitive to understand as the superquadrics.In the opinion of the domain experts, the reason for this lies, however, only in the fact that the local change of curvature is more sophisticated information with a broader space of possible variations.Therefore, it might even be interesting to modify the design of these glyphs to adapt for context or application-specific usage.

Conclusion and future work
This work presents a visualization approach for the comprehensive analysis of interface deformation in multiphase flows.The metric tensor contains information on the directions of stretching and their corresponding magnitudes, providing valuable insights into interface instabilities.Since the metric tensor does not capture the actual shape change, we have introduced a technique for the analysis of the curvature change that shows the bending strength and direction and hence allows for a comprehensive investigation of multiphase flows in terms of interface deformation.For visualization, we introduced a novel glyph design for stretching and bending.We then applied our method to visually analyze simulations with and without considering Marangoni convection, providing a useful tool to the domain experts for analysis.
In the future, we would like to apply our method to other scenarios and more complex simulations.Our visual analysis approach can support domain scientists who conduct state-of-the-art research on surface tension force, e.g., researching the influence of temperature differences, electromagnetic fields, and surface active agents (surfactants) on interface deformation.We would also like to investigate the applicability of the curvature change analysis to general mesh deformations.For instance, the computed curvature change (i.e., direction and strength) could be used to constrain mesh deformation.

Fig. 1
Fig.1Parabola-fitting at the interface for an example 2D case.a PLIC planes and their barycenters in a 3 9 3 neighborhood with resulting PLIC space for the central cell.b Rotated PLIC space and c parabola fitted to the interface barycenters

Fig. 2
Fig. 2 Computation of the curvature change at the interface.a Selected interface fragment (orange box).b Displacement vectors at the barycenters.c Resulting barycenter positions (red).Note that the velocities displayed are ũðx i Þ as calculated in Eq. 5. d, e Reconstructed paraboloids for the initial and advected barycenter positions, with corresponding principal directions k 1 and k 2 .f Computed paraboloid difference reveals the bending directions as the principal directions and their strength as the principal curvatures of the resulting paraboloid (with blue indicating bending downward and red upward)

Fig. 3 Fig. 4
Fig. 3 Glyph for interface deformation.a Glyph consists of two separate glyphs for stretching (red) and bending (green).Stretching is visualized using a disc-shaped glyph scaled along the stretching directions by their respective magnitude.Bending glyphs represent polynomials that indicate a curvature change in the principal curvatures' directions.b Applied to a droplet, color and opacity highlight regions of strong deformation.c Zoom-in allows users to analyze the deformation process, focusing on glyph shape. g

Fig. 5
Fig.5Different time steps in the Colliding Drops dataset (top row, a-f), and the Raindrops dataset (bottom row, g-l).The two water droplets in the Colliding Drops dataset have opposing velocities and collide slightly off-center, forming a disc that eventually disintegrates.Both drops in the Raindrops dataset move in the direction of gravity g at different speeds.The larger one is faster and eventually catches up with the smaller droplet, leading to a complex collision event and disintegration

Fig. 7 Fig. 8 a
Fig. 7 Stretching glyphs showing vertical stretching and horizontal contraction in the Raindrops dataset for a time step close to Fig. 5j.This stretching and contraction lead to the cohesion of the lower part, which can be observed in the subsequent time steps below

Fig. 10
Fig.10Comparison of the simulation where convection (I),d was considered in the solver, and one (II),e where it was disabled.In the case of the simulation without Marangoni convection, the neck between the two droplets grows almost symmetrically.This is reflected in both stretching and bending: in the vertical direction stretching, and a decrease in convexity can be observed.In the horizontal direction, there is an increase in convexity in the center, and to the sides at the base of the connection, there is an increase in concavity.In the case of (I),d Marangoni convection, there is less growth but an additional region on the left droplet exhibiting strong deformation.The white arrow in d indicates the direction of this capillary wave, as also indicated for a zoom-in on (I) c in Fig.11d

Fig. 11
Fig. 11 Close-up of the capillary wave propagating from right to left onto the left droplet (see blue arrow in d), with deformation glyphs showing interface stretching and bending.The bending glyphs show white strips to indicate the major directions (cf.Fig 4b).The stretching glyphs are extended with colored strips to highlight major directions and magnitudes (cf.Fig. 4c; color from contraction to stretching).Depending on simulation time, up to four different regions around the capillary wave can be observed in d: À in front of the wave, `on the rising slope, ´on the falling slope, and ˆbehind