{'fraction_non_alphanumeric': 0.025811442838555403, 'fraction_numerical': 0.006910011604599641, 'mean_word_length': 4.810092961487384, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An Adaptive Variational Fully Eulerian Fluid-Structure Interaction via Allen-Cahn Phase-field Modeling

In this study, we present a novel interface-driven adaptive technique using a fully Eulerian description for fluidstructure interaction (FSI). The proposed fully-Eulerian approach involves a fixed background unstructured mesh on which the fluid-structure interface is treated implicitly. The fluid-structure coupling is modeled by the phasefield finite element formulation based on the integration of the convective Allen-Cahn equation with the unified momentum equation for both solid and fluid dynamics. A bounded and stable solution of the Allen-Cahn equation is ensured via the positivity preserving variational formulation. To evaluate the solid stresses, the left CauchyGreen deformation tensor is convected at each time step to trace the evolution of the solid strain in the Eulerian reference frame. The adaptive refinement/coarsening for the unstructured mesh is determined by appropriate indicators and carried out by the newest vertex bisection algorithm. The proposed nonlinear adaptive partitioned procedure aims at reducing the amount of coarsening while simultaneously ensuring convergence properties of the nonlinear governing equations. We perform a detailed convergence and accuracy analysis via two benchmark problems namely, the pure solid system and a coupled fluid-solid system with an interface in a rectangle domain. We next systematically assess the performance of the adaptive procedure for the increasing complexity of problems. Finally, we demonstrate our fully-Eulerian interface-driven adaptive FSI model to simulate the contact problem between an elastic solid and a rigid wall using the well-known bouncing elastic ball problem.


Introduction
Fluid-structure interaction (FSI) is one of the most omnipresent and complicated problems in nature and engineering applications. Two-way interactions between a flowing fluid and a deforming solid occur in applications ranging from industry such as the occurrence of flutter in an airplane wing and flowing fluids in pipelines to biomedical fields such as flow of blood in arteries and bio-locomotion studies of fish, jellyfish, bat and among others [1,2]. The relative motion between a solid and the fluid in both internal and external flows leads to the application of loads on the solid resulting in deformations and/or displacements. These structural changes affect the fluid motion resulting in a two-way coupling between the solid and the fluid. This coupling is inherently nonlinear in nature, further complicating the dynamics of the problem for analytical methods [2]. The motivation for this work primarily comes from bio-inspired locomotion. For example, octopus-based soft robotics [3,4] involve very flexible structures (e.g., muscular hydrostats) with large deformations and self-contact which pose serious challenges to the mesh-based moving boundary and Eulerian-Lagrangian methods [2]. There is a need for an alternative approach to deal with these specific types of FSI problems that involve fairly complicated structural motions and constitutive relations [5]. In the present work, we explore the fully Eulerian approach to simulate such FSI problems with reasonable accuracy of the coupled dynamics.
For large-scale numerical FSI modeling, the continuum hypothesis for the physical domains is generally adopted. Using the continuum approach, the key challenges in modeling fluid-structure interactions pertain to the treatment of the boundary conditions at the fluid-solid interface; handling of dissimilar coordinate frames for the fluid and solid domains; and the construction of efficient and robust numerical techniques for the coupled nonlinear PDEs. While the mesh moves with the simulated continuum in Lagrangian description, the physical field/ continuum moves through a stationary or fixed grid in Eulerian description. Based on the choice of coordinate frames, one can consider fully Lagrangian, fully Eulerian, or hybrid Lagrangian-Eulerian approaches for solving fluid-structure interaction problems. In a Lagrangian approach, the fluid-solid interface is defined explicitly on a boundary conforming mesh following the motion of the solid in the Lagrangian coordinates, such as fully Lagrangian [6] and Arbitrary Lagrangian-Eulerian (ALE) [7,8,9]. Body-fitted ALE techniques have long dominated the realm of numerical solutions for FSI problems. While these techniques have been very successful for accurate interface tracking of several complicated multiphase and multiphysics engineering applications, they face difficulties dealing with flexible multibody contact dynamics (e.g., self-contact between octopus arms) or those involving large topological changes (e.g., the rupture of an aortic aneurysm).
Hybrid Lagrangian-Eulerian and fully Eulerian approaches attempt to address these difficulties. These techniques rely on a non-conforming discretization where the interface embeds or cuts across the fixed background grid. Such treatment allows handling large deformation problems in a relatively simple manner by removing the need for frequent remeshing required in body-conforming methods. Hybrid techniques such as the immersed boundary method (IBM) by Peskin [10,11] make use of the regularized delta functions to describe the interface, hence smoothing out stress discontinuities that are usually associated with fluid-solid interfaces. In immersed interface methods (IIM) [12], the regularized kernels are avoided and the jump conditions are directly built into the finite difference approximations. IIMs provide higher accuracy, although requiring coupling conditions that involve complex geometrical computations. Fedkiw et al. [13] worked extensively on the ghost fluid method, which is another example of a hybrid Lagrangian-Eulerian technique. Boffi et al. [14] and Zhang et al. [15] improved upon the traditional finite difference-based immersed boundary method by applying stabilized finite element technique to solve for both the fluid and the solid domains.
Some hybrid approaches introduce additional interface variables and couple the fluid and solid system using Lagrange multipliers e.g., fictitious domain method. In the fictitious domain approach, the fluid variables are extended inside the solid domain as well and the coupling between the two is carried out using a Lagrange multiplier. This multiplier imposes the kinematic constraint and in itself represents the body force to enforce this constraint. Glowinski et al. [16] used the fictitious domain approach to perform numerical simulations of fluid flow over moving rigid bodies. One can generate the solid mesh independently from the fluid mesh and the interface cuts the elements without the need for mesh alignment. There are variants of the fictitious domain approach such as the finite cell method which utilizes the high-order finite element and the framework of cut cells [17,18]. The challenging aspects of the finite cell method are the robust treatment of complex cut elements and the numerical integration over the cut cells. Similar to the fictitious domain method, the Eulerian-Lagrangian coupling methods based on the eXtended Finite Element Method (XFEM) do not require conformity between the meshes. The XFEM employs enriched basis functions that are designed according to the known solution form to improve the quality of the numerical solution. Chessa and Belytchsko [19] applied the XFEM to solve two-phase immiscible flow problems. Using the XFEM, one can simplify the interface discontinuity without conforming the mesh to the interface, hence removing the need for remeshing. Another method analogous to the XFEM, is the generalized FEM or GFEM based on the Galerkin method and partition of unity [20,21]. Similar to the finite cell method, the key problem associated with XFEM is the evaluation of integrals over discontinuous integrands on polyhedral domains. Vaughan et al. [22] carried out a comparison of the extended finite element method with the immersed interface method for elliptic equations. Felippa et al. [23] presented a classification of several interface coupling-based approaches for solving FSI problems. The majority of issues of the hybrid techniques arise from the coupling conditions between the Eulerian and Lagrangian variables. Methods such as immersed boundary, immersed interface as well as XFEM may suffer from inaccuracies for large topology changes and contact problems. In particular, some of these methods pose challenges in handling geometric complexities and numerical implementation. Techniques such as the finite cell method or XFEM can face the issue of extremely small cut cells which cause a problem during the linear solver step, thus requiring further rectifications [24]. A fully Eulerian approach is usually desirable for problems involving very large deformations of the solid.
As the name suggests, in a fully Eulerian approach both the solid and the fluid are represented in the Eulerian coordinates. This approach belongs to the fixed-mesh interface-capturing approaches. However, solid mechanics is formulated in terms of Eulerian coordinates in contrast to immersed and fictitious domain methods. The coupled system can be computed on a single mesh and the coupling can be constructed in a consistent variational monolithic formulation. Liu and Walkington [25] attempted to develop the evolution equation for the strain in an Eulerian framework for the flowing fluid containing visco-hyperelastic particles. The idea of a fully Eulerian approach for solving FSI problems was then demonstrated by Dunne [26], whereby the authors proposed an Eulerian framework for modeling FSI of an incompressible fluid and an elastic structure. The authors used a technique the so called the initial position set to track the initial location of a point. While this approach has been further extended by [27] for unsteady problems, Richter [28] used the present method for solving the non-linear interaction of an incompressible fluid and a hyperelastic solid. In [28], characteristic functions were introduced to identify the fluid and solid domains and the inverse map function is employed to trace the material coordinates at each time step. The inverse map function along with the characteristic function serves the purpose of the initial point set technique used in previous studies. In both works [27] and [28], the authors suggested the need for adaptive mesh refinement to maintain convergence in this approach. During the solution procedure, the evolution of the solid strains can be performed by convecting the inverse map function [29,30], the deformation gradient tensor or the left Cauchy-Green deformation tensor [25,31].
Handling the fluid-solid interface is one of the primary challenges in solving FSI problems. Hence, a substantial amount of effort has been devoted to perfecting this aspect of the numerical model by the FSI community. There are two key issues involved, namely: (i) accurate interface description i.e., satisfying the boundary conditions along the interface, and (ii) evolution of the interface in space and time. One way to capture the interface is via the sharp-interface approach (e.g., level set technique) where the fluid-solid interface is treated as a sharp boundary between the fluid and solid domains. The level set method was introduced by Sethian [32] and Fedkiw et al. [33] for tracking physical interfaces in Eulerian coordinates. The level set function is usually a signed distance function, which is used to describe the physical domains. Sussman et al. [34] introduced a reinitialization equation to retain the signed-distance property of the level-set function as it evolves in time. The reinitialization process may again prove to be computationally expensive and lead to poor mass conservation properties. Higher order discretisation methods and improved re-initialization [35,36,37] introduced for correcting the mass conservation issue for the level-set technique further make the process complicated and computationally expensive.
While the interface is defined sharply in the level set method, during its evolution in the Eulerian grid, the interface is subjected to numerical diffusion. This develops an uncontrollable diffused interface in actual simulations. Whereas in the diffuse-interface representation, we allow for diffusion in the model itself with the phase indicator (φ) varying smoothly across the interface (as shown in Fig. 1). Despite the lower accuracy in some cases, this approach provides more robustness and numerical stability. This is especially useful in problems involving large topological changes where interfacial dynamics play a significant role. The idea of diffused interfaces was introduced by van der Waals [38]. This was later taken forward and used for solving phase separation problems, where the Cahn-Hilliard [39] and Allen-Cahn [40] equations were established. The use of phase-field methods emerged in the 1990s for numerical modeling when they were applied to solidification processes [41,42]. We use the convective form of the Allen-Cahn equation to evolve the interface in space and time, in conjunction with the Lagrange multiplier term for mass conservation and positivity preserving terms [43] for a stable and bounded solution. The method is free from the reinitialization process from the level set method. Its free energy minimization process ensures energy stability and makes it more compatible with variational discretization. As a result, we prefer the phase-field method for interface capturing.
The present study builds upon our previous work [44], where an adaptive refinement procedure was presented for two-phase flows, integrating the Navier-Stokes and Allen-Cahn equations. The recently proposed positivity preserving variational (PPV) technique [43] is employed to maintain the stable and bounded solution of the coupled nonlinear differential system. The present work focuses on the successful integration of the unified continuum equations for the entire domain and the Allen-Cahn equation for the phase-field function and accurate implementation of the adaptive procedure for capturing the fluid-solid interface. Since we do not explicitly solve for the solid displacements, we evolve the solid strains to evaluate the deformations and stresses in the structure accurately. An adaptive algorithm based on the newest vertex bisection method [45] is employed for unstructured triangular meshes. This algorithm avoids the data structure to store the information about the refinement and coarsening nodes. The convergence criteria reduce the amount of coarsening while maintaining the convergence properties of the coupled discretized equations for the fluid-structure interaction. This article focuses on developing a fully Eulerian approach for solving FSI problems involving large solid motion and deformation with contact dynamics. While the solid dynamics is represented in the Eulerian coordinates with convection effect, one can expect a reduction in accuracy and conservation properties due to numerical diffusion. Hence, we introduce interface-driven adaptivity to improve upon these attributes of the problem. We work with residual-based adaptivity for the Allen-Cahn equation to specifically refine around the fluid-solid interface, hence reducing our computational efforts. To the best of the authors' knowledge, adaptive mesh refinement on an unstructured grid for fully Eulerian-based FSI applications is not carried out in the literature. We examine the defining parameters of the problem such as interface thickness parameter (ε) for different initial grid size h init and error tolerances through the use of performance indicators such as mass loss error (e mass ), elapsed time for the simulation and the residual error indicator (η). The primary focus of the work is to study the use of the diffuse interface approach for solving fluid-structure interactions via the fully Eulerian method and evaluate the adaptive refinement procedure for improved accuracy and conservation properties.
The article is organized as follows. Section 2 builds upon the mathematical formulation of the fully Eulerian approach by deriving the equations for the individual fields in the Eulerian frame of reference and further describes the associated equations. Section 3 discusses the semi-discrete variational forms for all the equations involved. Section 4 begins with a discussion on the need for adaptive refinement in a fully Eulerian approach. It further elaborates upon the procedures used to carry out mesh refinement adaptively. Section 5 presents a systematic convergence and validation study with the increasing complexity of problems. Finally, we apply the developed approach to a couple of application problems in Section 6 to demonstrate the efficacy and robustness of the solver. Section 7 summarises the key findings and discussions from the paper.

Fully Eulerian Model
This section deals with the development of the fluid and solid equations in the Eulerian frame finally culminating in the one-field formulation. From continuum mechanics conventions, let Ω denote the domain in material coordinates X at t = 0 and Ω t denote the deformed state of the domain in spatial coordinates, x. The domain is divided into two parts: the fluid and the solid, such that Ω = Ω f t ∪ Ω s t and ∅ = Ω f t ∩ Ω s t .

Fluid Equations in Eulerian Frame
The fluid domain is naturally expressed in the Eulerian coordinate system. The incompressible Navier-Stokes equations in the spatial domain Ω t can be written as: where the stress tensor is the rate of strain tensor. Velocity and pressure are a function of the spatial coordinates in the Eulerian frame, v f = v f (x, t) and p f = p f (x, t). When a fluid flow interacts with a structure, the fluid loading changes the configuration of the structure by inducing deformations and/or displacements. Structural response can be treated in Lagrangian or Eulerian coordinates.

Structure Equations in Lagrangian Frame
We assume the solid to be incompressible, hence all deformations are isochoric in nature. The governing equations for the structure in the Lagrangian coordinates are given bŷ where (·) represents quantities in the material configuration,P is the first Piola-Kirchoff stress tensor, and F is the deformation gradient tensor given by F = (I +∇u s ). Here u s denotes the displacement vector and I is the identity matrix.P is related to the Cauchy stress tensor (σ s ) by the relationP = Jσ s F −T and the material density is given byρ s = Jρ s , where J = detF is the volumetric stretch which is unity for incompressible materials, as mentioned above.

Structure Equations in Eulerian Frame
The above Lagrangian equations are converted to the Eulerian frame of reference using the deformation gradient tensor and the inverse map function (ξ(x, t) = X). We use the incompressible neo-Hookean model to describe the solid stresses in the present work. The conversion of the time derivative from Lagrangian to Eulerian coordinates can be obtained asρ where J is the volumetric stretch defined in the previous sub-section.

Remark 1.
It is worth noting that the conversion of these equations into the Eulerian frame introduces convective terms in the structural equation as well. This marks the primary point of difference compared to the Lagrangian approaches. These convective terms might cause issues during the numerical solution if not discretized appropriately with stabilization schemes.
Similarly, the stress term can be written as∇ Incompressibility condition of detF = 1 in the Lagrangian sense translates to ∇ · v s = 0 in the Eulerian frame [46]. Thus, the governing equations for the solid domain in the Eulerian frame are given by Material configuration Spatial configuration Figure 2: Representation of the inverse map function (ξ) that maps from the current (spatial) configuration to the reference (material) configuration where the Cauchy stress tensor can be expressed as In the above equation, σ s sh = µ s L (F F T − I) in accordance with the incompressible neo-Hookean model (refer to Appendix B for discretisation of the solid shear stresses). In the previous relation, µ s L is the shear modulus of the material (Lame's second parameter), given by µ s L = E 2(1+ν) , where E is the Young's modulus of the material and ν is the Poisson's ratio.

Remark 2.
The presence of the viscous term in Cauchy stress tensor equation allows us to model visco-hyperelastic solids as well or simply to stabilise the numerical scheme in the absence of visco-elastic solids (it can be set to zero for all other cases).

Evolution of the Solid Strain
Apart from solving the momentum equation for the solid, it is also essential to evolve the solid strain to evaluate the deformations and stresses correctly. This becomes crucial since we are not directly computing the solid displacements. This additional information closes the system of equations for the fully Eulerian model. The mapping from the spatial to the material coordinates is carried out through the inverse map function ξ such that X = ξ(x, t). This mapping function helps to trace the material points on the solid during its motion, as shown in Fig. 2.
Consider initial undeformed configuration of the solid ξ(x, 0) = x = X. By convecting the inverse map function [29,30], the evolution of ξ subject to this initial condition can be written as: In the continuum sense, one can show the equivalence between the the convection of inverse map function and the deformation gradient tensor or left Cauchy-Green deformation tensor [25,31] using simple manipulations and the fact that F = (∇ξ) −1 and B = F F T (refer to Appendix A). The transport equations for the deformation gradient tensor and left Cauchy-Green tensor can be expressed as follows: and where F is the deformation gradient tensor and B is the left Cauchy-Green deformation tensor. The left Cauchy-Green tensor equation is chosen for the evolution of the solid strains for two important reasons: • If we solve for ξ, we would still need to compute either F or B to evaluate the solid stresses. Since these quantities are not directly available, we will be required to calculate the gradient of ξ and use some interpolation/projection technique to obtain the nodal values at every time step. These additional computations can be time-consuming and may introduce new errors.
• The preference of B over F arises from the simple fact that B is a symmetric tensor.
Due to the application of the unified approach, there is a need for an additional variable to keep track of the fluid-solid interface, which is described in the following section.

Interface Representation via Phase-field Method
In this section, we focus on an important aspect of any FSI problem i.e. dealing with the fluid-solid interface. As mentioned previously, we use the diffused interface based phase-field method for capturing the interface. In the phase-field representation, the continuum is considered as a mixture of the two phases and free energy is stored. Bulk phases are regions where there is a slow variation of the phase-field function and interface is the region where there is a high localized variation of the phase-field function. This avoids the requirement of the knowledge of exact interface location at every time step, further preventing setting boundary conditions at the interface for the bulk phases. This is one of the primary advantages of the phase-field method over sharp-interface methods such as the level set technique. In both Cahn-Hilliard and Allen-Cahn equations for phase-field description, the phase separation and phase transition in the diffused interface result from the process of free energy minimization, where the Ginzburg-Landau free energy functional is given by: where Ω is the bounded physical domain, H 1 (Ω) denotes the space of square-integrable real-valued functions with square-integrable derivatives on Ω, L 4 (Ω) denotes the function space in which the fourth power of the function is integrable, R ≥0 represents the set of non-negative real numbers and φ(x, t) is referred to as the order parameter or the phase-field function which indicates the two domains and has values of either +1 or -1. The first term in the RHS is the bulk or mixing energy and depends on the local composition of the mixture. The second term is the interfacial or gradient energy and depends on the composition of immediate environment. Hence there is a trade-off between preference of pure phases and a mixed uniform phase. The ratio of these two effects controlled by ε decides the thickness of the diffused interface region. At equilibrium, the interface thickness is the distance over which φ varies from −0.9 to 0.9 which can be estimated as 2 The Allen-Cahn equation is computationally preferred for its lower order equations compared to the Cahn-Hilliard equation. The missing mass conservation property in the Allen-Cahn equation is made up for by adding a Lagrange multiplier [47,48] or an anti-curvature term [49]. The final convective form of the Allen-Cahn equation after addition of the Lagrange multiplier can be derived as follows: where β(t) is the time-dependent part of the Lagrange multiplier, given by . The functional form for the bulk energy F (φ) has been chosen as the double-well potential function F (φ) = 1 4 (φ 2 − 1) 2 in the present study.
Remark 3. This quartic polynomial function is, however, an approximation to the more accurate logarithmic function that is derived from thermodynamic considerations. Analogous studies have been carried out for the more rigorous Cahn-Hilliard equation [50,51] to illustrate its numerical properties under the assumption of logarithmic free energy. The polynomial approximation, though has been shown to perform reasonably well in literature [52] and is more popular due to its numerical simplicity.
The boundary condition and initial condition for the Allen-Cahn equation are as follows: This completes the description of the mathematical formulation of the system. In the next section, we present the weak variational forms of the above equations.

Semi-discrete Variational Eulerian FSI Formulation
In the earlier section, we presented the mathematical models and the corresponding governing equations representing them. Now we will move on to discretizing those non-linear coupled equations in order to be able to solve them numerically. This section presents the semi-discrete variational forms of the governing equations using the finite element method with additional stabilization terms wherever required. The temporal discretization is carried out using the generalised-α method [2,43].

Unified Continuum Equations
To begin with, we present the unified continuum equations for the fluid-solid system. Combining the developments from the previous sections, we arrive at the one-field continuum equations for the fluid and solid domains using the phase indicator φ(x, t) as where the phase properties are defined as follows: Hence the fully Eulerian formulation results in a very simplified and robust set of equations for the combined fluid-solid domain. Let S h be the space of trial solutions, whose values satisfy the Dirichlet boundary conditions and V h be the space of test functions whose values vanish on the Dirichlet boundary. The variational form for the unified continuum equation can be written as: The first line consists of the Galerkin terms for the combined momentum and continuity equations. The second line contains the Petrov-Galerkin stabilization terms for the continuum equations. R m and R c are the elementwise residuals for the momentum and continuity equations respectively. The stabilization parameters τ m and τ c [53,54,55] are given by where C I is a constant derived from the element-wise inverse estimate and G is the element contravariant metric tensor.

Left Cauchy-Green Tensor Transport Equation
We next proceed to the variational form of the B tensor equation to evaluate the solid stresses accurately. Let S h be the space of trial solutions, whose values satisfy the Dirichlet boundary conditions and V h be the space of test functions whose values vanish on the Dirichlet boundary. The variational form of the left Cauchy-Green tensor equation can be written as: where ':' represents component wise multiplication. We obtain updated values of B by solving the above variational equation which can then be used to calculate the new solid stress values before substituting in the unified continuum equations in the next time step.

Allen-Cahn Equation
Finally, we present the interface modeling aspect of the formulation. The semi-discrete Allen-Cahn equation in the convective-diffusive-reactive form can be written as: where Convection velocity = v, The above expressions have been obtained by simplifying terms using the generalized-α time integration [43]. Let S h be the space of trial solutions, whose values satisfy the Dirichlet boundary conditions and V h be the space of test functions whose values vanish on the Dirichlet boundary. The final stabilized form of the Allen-Cahn equation can be stated as: The first line of Eq. (19) contains the Galerkin terms, the second line contains the linear stabilization terms with the stabilization parameter, τ [53,55] given by The third and fourth line consist of the positivity preserving terms. The positivity condition is enforced to prevent oscillations in φ at regions of high gradient which might have otherwise led to negative values of density or viscosity. Since upwinding makes the solution first order accurate, hence the addition of diffusion is restricted to regions of oscillations only. To achieve this criteria, the non-linear positivity preserving terms are made to depend on the residual of the equation [43]. The implementation of the PPV method makes the scheme at least second order accurate in space and the method is able to capture the high gradient internal and boundary layers in multi-dimensions. The χ |R(φ h )| |∇φ h | factor adds a nonlinear property and provides the functionality of a limiter for the upwinding near oscillating regions. Now we proceed to the central topic of this paper i.e. adaptive mesh refinement. We present the adaptive refinement technique in the next section, based on the aforementioned developments of unified continuum equations and positivity preserving Allen-Cahn formulation.

Adaptive Variational Formulation
In this section, we discuss the adaptive mesh refinement technique for the coupled solver. Body-fitted techniques such as ALE solve for the solid displacements exactly and track the interface in an explicit manner. Lack of exact solid displacements in an Eulerian framework and use of a diffused interface approach tend to bring down the accuracy of the solution. Hence it becomes essential to use local refinement techniques for complicated problems for preservation of accuracy and conservation properties. Since the fluid-solid interface is the source of maximum error in the solution, it is naturally wiser to use interface driven adaptivity techniques for the mesh to improve overall convergence and accuracy of the numerical procedure and reduce computational effort.
In this work, this has been carried out via residual dependent adaptivity which is based on the reasoning that at locations where the solution exhibits oscillations, the value of the residual is usually higher. This occurs due to dominating convection and reaction effects on an under-resolved mesh. Hence our goal is to minimize this residual to obtain accurate and converged solution. With that goal in mind, we first define some of the errors that characterize the current numerical framework. The velocity, pressure and order parameter increments are evaluated in each time step using Newton-Raphson type iterations. The linearized system for the unified continuum equations can be expressed by: where K Ω is the stiffness matrix of the momentum equation consisting of inertia, convection, diffusion and stabilization terms, G Ω is the discrete gradient operator, G T Ω is the divergence operator and C Ω is the pressure-pressure stabilization term. Here, ∆v and ∆p are the increments in velocity and pressure, respectively and R m (v, p) and R c (v) represent the weighted residuals of the stabilized momentum and continuity equations respectively. Let the updated quantities at t n+1 (k) be represented as X n+1 (k) , k being the non-linear iteration index and the increments in these quantities be represented as ∆X. The error in solving the unified continuum equations can be defined as Similarly, the linearized form for the transport of the left Cauchy-Green tensor can be written as: where K CGT is the stiffness matrix and R(B) represents the weighted residual for the left Cauchy Green tensor equation. The numerical error in solving the above equation can be written as Finally, the linearized form for the Allen-Cahn equation can be expressed as: where K AC consists of the inertia, convection, diffusion, reaction and stabilization terms and R(φ) represents the weighted residual for the stabilized conservative Allen-Cahn equation. Similar to the error evaluation for the deformation tensor, the numerical error in solving the Allen-Cahn equation can be written as Based on discussions in Joshi and Jaiman [44], the adaptive error indicator can be derived from the discretized Allen-Cahn equation. Consider the domain Ω which consists of elements Ω e , chosen such that Ω = ∪ n el e=1 Ω e and ∅ = ∩ n el e=1 Ω e where n el is the number of elements. Let Γ be the Lipschitz continuous boundary of the domain Ω, Γ D and Γ N be the Dirichlet and Neumann boundaries of Ω respectively such that Γ = Γ D ∪ Γ N . Furthermore, let E denote the set of edges for all the elements in the domain, E Γ , E Γ D and E Γ N be the set of edges on the boundary Γ, Γ D and Γ N respectively. Consider E Ωe to be the set of edges of an element Ω e and E Ω be the set of all the interior edges of Ω. The Galerkin terms of the variational Allen-Cahn equation can be written after integration by parts as where R Ωe and R E are the element and edge based residuals given as, where J E (·) is the jump of the argument across the element edge E and n E is the normal to the edge E. After a few algebraic manipulations and using the Cauchy-Schwarz inequality [44], we can define the error estimate between the true (φ) and approximate (φ h ) solutions using the following relation: where In the above equations, η is the error estimate, η Ωe is the error indicator for a single triangular element Ω e , h Ωe is the diameter of the element Ω e and h E is the length of the edge E. This is a better approach, as here we do not need to solve any additional equation i.e., the 'dual problem' [30] and can directly obtain the error residuals from the solution of the Allen-Cahn equation.
After obtaining the error estimates from the Allen-Cahn equation, the newest vertex bisection algorithm [45] is employed to carry out the adaptivity procedure. The algorithm is based on appropriately marking the elements and edges that need to be coarsened or refined and then carrying out the required procedure. The Dörfler criterion [56] is used to mark the elements for refinement which finds the minimum set of elements Ω M such that where θ is a user defined parameter in the range (0, 1). This criterion ensures that the elements which have a larger contribution to the error estimate are chosen for refinement.
The newest vertex bisection algorithm is used for refining the triangular elements of the unstructured mesh along the diffused interface as shown in Fig. 3. The longest edge (BC) of the element ABC is chosen first and Figure 3: Application of the newest vertex bisection method for refinement of a portion of the representative unstructured grid along the diffused interface: ∆ ABC has been marked for refinement by the Dorfler criterion. Nodes D, E and F are mid-points of the edges to which bisectors are drawn to create additional elements, starting with node D since BC is the longest edge. In the above figure, δ ≈ 4ε at equilibrium (where ε is the interface thickness parameter). bisected by placing a node (D) at the mid-point. Then the remaining edges are bisected by placing nodes E and F at their mid-points respectively, in accordance with the corresponding edge being opposite to the newest vertex created in the respective elements. Once the particular element is refined, the hanging nodes are avoided by the creation of additional elements adjoining the triangle. The same procedure is then repeated for all the marked elements along the interface. An additional array is created in the process which maintains the identities of the nodes added by bisection and those that belong to the initial grid. The nodes added later are eligible for removal during the coarsening procedure, while initial nodes cannot be coarsened. This step considerably simplifies the numerical implementation and helps in the coarsening step as well. Once all the marked elements are refined, the governing equations are solved on the refined mesh till convergence. Elements with very small residual errors can be coarsened using the coarsening algorithm. This is done to reduce the computational effort for the problem. Similar to the refinement criterion, a Dörfler criterion is defined for coarsening as well, with θ c as the user-defined coarsening parameter. The earlier created array that contains the information about the added nodes by bisection helps in identifying the nodes that are suitable for coarsening [57]. Hence we cannot coarsen below the initial number of elements that we start with. None of the coarsening/refinement steps are carried out in the intermediate iterations while implementing the algorithm. This ensures that the non-linearities of the governing differential equations are captured.

Algorithm
The proposed algorithm for carrying out adaptive procedure for a phase-field based fully Eulerian FSI solver has been presented in Algorithm 1. At a given time step, we start with all the primitive variables: velocities v n , pressure p n , order parameter φ n and left Cauchy-Green tensor B n at the nodal points on grid T n . Before entering the non-linear iteration loop, we predict the solutions for the next time step: v n+1 (0) , p n+1 (0) , φ n+1 (0) and B n+1 (0) . Inside the non-linear iteration loop, we first solve the unified continuum equations to obtain the velocities and pressures in the entire domain as well as evaluate the corresponding error e U C . The updated velocities are then used for evolving the solid strain by solving the transport equation for B. The corresponding error e CGT is evaluated at this step. Finally, we solve the Allen-Cahn equation to evolve the fluid-solid interface in space and time and estimate the error e AC .
Next, the error estimator for the adaptive procedure is computed using Eq. (30) and Eq. (31). The succeeding steps are used to decide if the present triangulation (mesh) needs refinement/coarsening based on the user-defined tolerances for the convergence criteria. nIterM ax represents the maximum number of nonlinear iterations and nElemM ax represents the maximum number of elements allowed in the domain. The boundary conditions are satisfied at the end of each loop (line 27) and the updated solution values are copied to the previous time step variables upon exiting the non-linear loop (line 31).
Algorithm 1 Adaptive procedure for a fully Eulerian FSI solver 1: Given v 0 , p 0 , φ 0 on T 0 2: Initial refinement step 3: Loop over time steps n = 0, 1, · · · 4: Start from known variables v n , p n , φ n , B n on T n

Implementation Details
The above algorithm is implemented on an unstructured mesh imported from Gmsh [58]. The CAD drawing and the meshing are carried out via Gmsh toolbox. The mesh file is first pre-processed to separate the nodal data, connectivity data, and boundary data. This step is followed by the initialization of velocities, pressure, order parameter and deformation tensors in the domain and the application of boundary conditions. The present triangulation is then passed through an iterative initial refinement step before entering the time loop. This is done to start with a sufficiently refined mesh along the interface to avoid high gradients for velocities and stress tensors when we enter the time loop despite the coarse initial mesh. The user-defined parameters used for this step are θ 0 = 0.4 and θ c,0 = 0.3. The system then enters the time loop where the nonlinear set of equations are solved in a partitioned manner [43]. Implicit discretizations are carried out for the underlying equations for improved efficiency and robustness of the solver. The incremental quantities are computed by employing Newton Raphson type iterations and the generalized-α time integration is used to update the solver. The Jacobians are calculated and the matrices are assembled for each of the equations separately and then solved sequentially as outlined in the previous sub-section. The updated velocities from the solution of the unified continuum equations are used for evolving the solid strains via the B transport equation and updating the order parameter values through the Allen-Cahn equation. The updated order parameter values are further utilized in the interpolation of density ρ(φ), viscosity µ(φ) and stresses σ(φ). These updated quantities and the solid stresses computed from B are fed into the unified continuum equation for the next iteration. Gaussian quadrature technique is used for the integration of the variational equations. The present code is vectorized to reduce the excess time spent in loops.
The solver is designed such that refinement of the mesh is carried out till the criteria for tol R is satisfied and coarsening is only carried out in the last non-linear iteration. The governing equations are then solved on the refined mesh to obtain updated solution values. This is done to reduce the individual errors of the corresponding equations below their tolerance values on the same background triangulation to capture the non-linearities accurately and ensure proper convergence properties of the governing equations. We also coarsen the grid if the number of elements n el exceeds nElemM ax. Based on the proposed procedure, we present several numerical tests of increasing complexity to demonstrate the effectiveness of the solver.

Convergence Studies
This section elaborates upon the convergence and error analysis of the fully Eulerian formulation. Two different test cases have been solved using the present solver with increasing order of complexity, namely a pure solid system and a fluid-solid interaction system. In each case, the convergence of the solver with respect to the Eulerian grid spacing is analysed. The section concludes with a popular validation case for demonstrating the correctness of the solver.

Pure Solid Wall
The first attempted benchmark problem is that of a pure solid system which was carried out to verify the structural components of the solver. The simulation is performed on a square domain consisting entirely of an elastic Neo-Hookean solid. The method of manufactured solutions is used to validate the solver against the reference velocity field whose streamfunction representation is as follows: ψ(x, t) = ψ 0 sin(ωt)sin(k x x)sin(k y y), where ψ 0 = 0.25, ω = 2π and k x = k y = 2π.
The external body force that is applied to reproduce the above motion is: The problem is solved on the above domain on various grids ranging from h = 0.083 to h = 0.01. The L2 norm of the difference between the expected (reference) and obtained velocities for different grid sizes is plotted in Fig.  4b. We observe second order convergence for this problem with spatial order of accuracy O(h p ) with p ≈ 2.2.

Flow Induced Deformation of an Elastic Wall
The second problem to be considered is the flow induced deformation of an elastic wall. This benchmark was carried out to verify the full FSI solver using the fully Eulerian approach. This is a modified lid-driven cavity problem where the bottom part is occupied by a soft elastic solid to capture the deformations in the domain. The stiffness of the solid is reduced greatly and convection is removed from the system to visualize the deformations accurately. The size of the domain is x ∈ [0, 0.3], 1, x ∈ (0.3, 1.7), sin 2 (π(x − 2)/0.6) x ∈ [1.7, 2]. The no-slip condition is imposed on the rest of the boundaries and Neumann condition is applied for the orderparameter on all the boundaries. The phase-field method has been used for capturing the interface on a structured grid first and later used on an unstructured grid with adaptivity. The initial condition for the phase-field function φ is given by, where ε is the interface thickness parameter. Figure 5 illustrates the order of accuracy plot and a comparison of the present results with those of [26] for the non-adaptive case. The accurate interface capturing for this benchmark problem ensures the applicability of the present method. A systematic convergence study is carried out for ε = 0.01 for different mesh sizes ranging from h = 0.16 to h = 0.005. The x-velocity solution on the finest mesh (h = 0.005) is taken as the reference solution. The L2 norm of the error between the velocity field for the given mesh with respect to the reference solution is plotted in Fig. 5c. The spatial order of accuracy comes out to be around O(h p ) with p ≈ 1.6. The decrease in the order of accuracy from previous case can be attributed to the presence of the fluid-solid interface.

Adaptive Investigations
We carry out the above problem with adaptive mesh refinement as well. Figure 6 illustrates the refinement of the mesh along the fluid-solid interface almost exclusively. A comparison between the adaptive and non-adaptive cases has been carried out in Fig. 7. It can be observed that the total mass loss error reduces by nearly 4.5 times for the adaptive case compared to the non-adaptive case. The mass loss error of the order parameter is evaluated as: where m is defined as m = Ω φ n dΩ and m t=0 is the mass of the order parameter at t = 0. The increased density of the elements helps in accurate capturing of the fluid-solid interface, hence resulting in reduced mass loss error. Similarly the cumulative elapsed time decreases by a factor of nearly three for the adaptive case thus resulting in reduced computational costs.

Effect of Interfacial Resolution
A parametric study was carried out to investigate the impact of interfacial resolution using the performance indicators of RMS errors for e mass and η (adaptive error indicator). The RMS value of any generic quantity λ  Figure 8 illustrates the variation of the RMS mass loss error with the interface thickness parameter ε. At very low values of the interface thickness parameter (Region I), we do not have enough grid cells to capture the interface. Even after adaptive refinement, we are not able to accurately describe the interface with such low ε values leading to high mass loss errors. On the other hand, when the interfacial thickness parameter is much higher (Region III), the interface is way too diffused. The grid refinement does not have a significant effect on the final resolution of the interface as we are constrained by the high ε values, again leading to high mass loss errors. Hence, from these plots we are able to come up with an optimum range for ε (Region II) in which simulations can be carried out. Based on Fig. 8, we can conclude that ε should lie between 0.004 − 0.014 to minimise the total mass loss error. The small value of this parameter can be attributed to the initial coarse mesh which is possible due to the adaptive refinement in the later time steps. Hence despite starting with a fairly coarse mesh, we recover similar levels of accuracy by specifically refining the grid along the fluid-solid interface. The RMS adaptive error indicator, η RM SE doesn't vary significantly with the given parameter (Fig. 9). Hence, we can conclude that the adaptive φ :

Soft Solid in a Driven Cavity
We simulate a soft solid undergoing deformation and motion inside a lid-driven cavity to further validate the solver. The size of the domain is [0, 1] × [0, 1] and the radius of the circular solid is 0.2 units, placed at (0.6, 0.5) initially. The physical parameters of the problem are: fluid viscosity µ f = 10 −2 , solid viscosity µ s = 10 −2 , shear modulus µ s L = 0.05, fluid density ρ f = 1 and solid density ρ s = 1. This problem has been solved previously using the finite difference [59] and finite volume [46] methods via the fully Eulerian approach. The initial condition for the order parameter is assumed to be, The solid is initially at rest inside the square domain and the top lid is given a horizontal velocity of 1 unit, V top = 1. No-slip boundary condition is imposed on the other three walls for velocity and Neumann boundary condition is imposed on all the boundaries for the order parameter. The solid moves due to the motion of the surrounding fluid and deforms accordingly. The hydrodynamic repulsion due to the presence of the liquid layer is able to prevent any unwanted overlap/penetration between the soft solid and the rigid wall. Figure 11 illustrates the location and deformation of the elastic solid inside the lid-driven cavity at various time instants. The black solid lines indicate the streamlines inside the domain. We have obtained a comparison of the trajectory of the centre of mass of the solid for various initial grid sizes with the reference solution of [59] in Fig. 12a. Finally, Fig. 12b illustrates the mass loss error of the system over time for the different cases considered.

Application Problems
The above-developed formulation is now applied to a more complicated problem of a ball dropping in a viscous fluid, involving contact to evaluate the robustness of the solver. The second example is an extension of this problem where two elastic balls are dropped from different initial heights. Elastic solid  Elastic ball Since, gravity is the only external force acting, the ball experiences "free-fall" inside the viscous fluid. The ball begins to move down due to the action of gravity with increasing speed. It comes in contact with the floor at point A (Fig.14b) and begins to compress. Due to the smooth, elastic nature of the ball, it stores its energy during compression and re-uses it to rise up during the rebound process. The important instances during the motion of the ball have been marked in Fig. 14b. The bottom-most part of the ball comes in contact with the floor for the first time at point A and compresses from A to B; the ball begins to rise at point B where we observe a change in direction of the velocity and finally, attains its maximum height at point C where the velocity again becomes zero.
We use adaptive mesh refinement for this problem as well (Fig. 15a). The interface thickness parameter ε was taken as 0.01, for an initial grid size of h init = 0.04 as suggested in the previous study. The user-defined parameters for the Dorfler criterion were taken as θ = 0.5 and θ c = 0.05. A significant improvement in the conservation properties is observed with the introduction of adaptive refinement in the grid as illustrated in Fig.  15b. This can be again attributed to accurate resolution of the fluid-solid interface with the introduction of the adaptive procedure.

Bouncing of Two Balls in a Liquid Tank
Next, we present a demonstration problem of multiple elastic bodies with free-motion solved using the phasefield fully Eulerian approach without adaptivity. Here we observe the motion of two elastic balls placed inside a liquid tank containing a viscous incompressible fluid. The radius of each ball is 0.3 units initially. The balls are released from different initial heights and are allowed to fall under the action of gravity. The initial condition for Similar to the previous case, the balls experience collision and rebound from the bottom wall by virtue of the stored elastic energy during the deformation process. The higher ball experiences stronger rebounce due to larger initial energy compared to the ball placed at a lower height. We also observe that the higher ball drifts away from the wall as it moves downwards. Figure 16b illustrates the contour plot for the magnitude of velocity at the final time step. Figure 17 demonstrates the motion history of the two individual balls. The right ball undergoes several minor rebounds from the floor of the tank, whereas the left ball undergoes one major rebound during the time of observation.
Adaptive refinement is carried out for this problem with ε = 0.01 for an initial grid size of h init = 0.05. The adaptive mesh at the start of the time loop is shown in Fig. 18a. Figure 18 illustrates the improvement in conservation properties with the introduction of adaptive mesh refinement. We observe significant reduction in error of individual ball volumes (Fig. 18b) and the total mass of the order parameter (Fig. 18c) for the adaptive case compared to the non-adaptive cases, owing to the increased density of elements at the fluid-solid interface. Elastic balls ||v|| : The presented implicit variational fully Eulerian algorithm is general in its formulation and can be extended to three dimensions, without any restriction from the proposed variational formulation. However, the geometric complexities involved in refining and coarsening of the grid in three dimensions pose some challenges during parallelization [60,61]. These are some of the aspects for future consider along with the incorporation of constitutive models for complex solid deformations and multibody contact problems.

Conclusions
In the current work, we presented an adaptive variational fully Eulerian approach for solving fluid-structure interaction problems on an unstructured grid to handle large-scale topological changes in the solid. Both fluid and solid domains are described on an Eulerian grid and the interface is represented using the phase-field formulation. The solid strains are evolved using the left Cauchy-Green deformation tensor to evaluate the solid stresses and de-formations. We introduced adaptive mesh refinement to improve accuracy and conservation properties originally lost due to the Eulerian representation of the solid. The adaptive error indicator relies on the residual of the Allen-Cahn equation and hence the procedure effectively results in an interface-driven adaptive grid. The refinement and coarsening steps are avoided in intermediate Newton iterations to capture the non-linearities of the governing differential equations. We presented several test cases with increasing order of complexity to demonstrate the effectiveness and robustness of the present approach. The second-order spatial accuracy of the solver was established through a systematic convergence study. The comparison studies yielded a nearly three times reduction in computational effort and about 4.5 times improvement in mass loss error for the adaptive cases compared to the non-adaptive ones. We performed a parametric study to suggest the optimum range for the interface resolution (ε) for different initial grid sizes to minimize the mass loss error. We demonstrated that irrespective of the initial grid size, we have a recommended range of the interface thickness parameter, ε = 0.004 − 0.014, to perform FSI simulations. We also performed a validation study of the present approach using the deformation of a soft solid inside a lid-driven cavity. Finally, we presented application problems to showcase the handling of contact and large-scale deformation in the final section.
In this appendix, we show that the evolution of the solid strains can be carried out by either convecting ξ, F or B and that the three methods are analytically identical to each other. As mentioned previously, the above transport equation requires that the material derivative of the inverse map function (ξ) is zero, subject to an initial condition: ∂ξ ∂t + (v s · ∇)ξ = 0, By taking the material derivative of the left cauchy-Green tensor B = F F T , we can write: Hence, the three equations of ξ, F or B for evolving the solid strains are equivalent in a continuum sense.

Appendix B. Discretization of Solid Shear Stresses
We expand the solid shear stresses in terms of the left-Cauchy Green tensor to be incorporated in the unified continuum equations. The shear stress in the solid, according to incompressible Neo-Hookean model, is given by: = B n + α∆t(∂ t B n + γ(∂ t B n+1 − ∂ t B n )), (B.2b)