An adaptive isogeometric analysis approach to elasto‐capillary fluid‐solid interaction

We present an adaptive isogeometric‐analysis approach to elasto‐capillary fluid‐solid interaction (FSI), based on a diffuse‐interface model for the binary fluid and an Arbitrary‐Lagrangian‐Eulerian formulation for the FSI problem. We consider approximations constructed from adaptive high‐regularity truncated hierarchical splines, as employed in the isogeometric analysis (IGA) paradigm. The considered adaptive strategy comprises a two‐level hierarchical a posteriori error estimate. The hierarchical a posteriori error estimate directs a support‐based refinement procedure. To attain robustness of the solution procedure for the aggregated binary‐fluid‐solid‐interaction problem, we apply a fully monolithic solution procedure and we introduce a continuation process in which the diffuse interface of the binary fluid is artificially enlarged on the coarsest levels of the adaptive‐refinement procedure. To assess the capability of the presented adaptive IGA method for elasto‐capillary FSI, we conduct numerical computations for a configuration pertaining to a sessile droplet on a soft solid substrate.


INTRODUCTION
Elasto-capillary fluid-solid interaction (FSI, or, concisely, elasto-capillarity) pertains to the deformation of an elastic solid due to capillary forces in a contiguous fluid meniscus. The fluid meniscus generally carries surface energy, which leads to a localized force on the solid surface at the contact line, and corresponding deformation of the elastic substrate. Elasto-capillarity is of fundamental importance in a variety of high-tech industrial applications, such as inkjet printing, additive manufacturing and lab-on-chip technologies, as well as in miscellaneous biophysical processes. Various complex physical phenomena emanate from elasto-capillarity, for example, durotaxis, 1 that is, seemingly spontaneous migration of liquid droplets on solid substrates with an elasticity gradient; capillary origami, 2 that is, large-scale solid deformations by capillary forces; the inverse cheerios effect, that is, long-distance attraction or repulsion of sessile droplets on deformable substrates; 3,4 and, in the case of visco-elastic solid substrates, stick-slip dynamics. 5,6 Elasto-capillarity moreover plays an important role in experimental investigations of the Shuttleworth effect, that is, deformation-dependent surface tension of solids. [7][8][9] Despite significant progress in computational models and simulation techniques for the interaction of solids and conventional (mono-)fluids (see References 10,11 for an overview) and for binary fluids separately (see, for instance, , modeling and simulation of elasto-capillary FSI has remained essentially unexplored. Recently, some first simulation approaches have appeared based on diffuse-interface models for the binary fluid; see  Diffuse-interface binary-fluid models yield a complete and versatile modeling paradigm for elasto-capillarity in view of their implicit treatment of topological changes of the fluid meniscus, which includes contact-line slip at the solid-fluid interface associated with dynamic wetting, and their rigorous thermodynamic basis. A comprehensive thermodynamic framework for elasto-capillary FSI based on a diffuse-interface binary-fluid model has been established in Reference 21. Diffuse-interface models have historically been propounded as a compelling alternative to conventional sharp-interface approaches for binary fluids, in view of their implicit treatment of the fluid-fluid meniscus. 23 In particular, in computational methods for binary-fluid flows, sharp-interface models require an explicit representation of the interface, for example, by a moving mesh that conforms to the meniscus or by means of a level-set representation, while diffuse-interface models provide an implicit representation of the meniscus by means of the order parameter. However, this perspective is one-sided: while, indeed, in diffuse-interface methods it is not necessary to track the meniscus explicitly, it is necessary to resolve it. In most applications, the thickness of the fluid meniscus is many orders of magnitude smaller than other relevant length scales in the problem under consideration, for example, the radius of a droplet. For instance, x-ray measurements 24 and molecular-dynamics simulations 25 have conveyed that the thickness of a liquid-vapor interface is typically a few Å, while the radius of a pico-liter droplet, which are the smallest droplets used in current technologies, is approximately 10 m. This amounts to a ratio of length scales of more than 10 4 , which is clearly beyond the resolution capabilities of conventional numerical approximation methods. For nano-and microliter droplets, the scale disparity is even more pronounced, increasing to 10 5 and 10 6 , respectively. Fortunately, for most purposes it is not necessary to adhere to the actual value of the interface thickness, and it suffices that the thickness of the diffuse-interface is sufficiently small (in some appropriate sense) relative to other length scales in the problem. For instance, elasto-capillary FSI of a sessile droplet with a soft solid substrate leads to the formation of a so-called wetting ridge in the vicinity of the contact line; see, for example, References 5,21. This wetting ridge culminates in a kink at the contact line. For diffuse-interface binary-fluid models, this kink in fact exhibits a finite radius of curvature, which vanishes in the sharp-interface limit, that is, as the thickness of the diffuse-interface passes to zero. 26 The overall deformation of the solid-fluid interface including the wetting ridge is however insensitive to the thickness of the diffuse interface, provided that the diffuse-interface thickness is sufficiently small relative to the ridge's characteristic dimensions. The dimensions of the wetting ridge are proportional (with a factor close to unity) to the elasto-capillary length scale composed of the ratio of surface tension to the shear modulus. Typically, the elasto-capillary length is on the order of 10 m, and it plausibly suffices to set the diffuse-interface thickness within the range of 10 −1 m and 1 m. For most applications, the ratio of the droplet radius to the thickness of this mollified diffuse-interface, however, still leads to prohibitive resolution requirements for nonadaptive approximation methods.
In this work, we consider an adaptive isogeometric analysis (IGA) approach to elasto-capillary fluid-solid-interaction (FSI) problems, based on an incompressible Navier-Stokes-Cahn-Hilliard (NSCH) diffuse-interface model for the binary fluid 15 and an Arbitrary-Lagrangian-Eulerian (ALE) formulation for the FSI problem. In particular, we consider approximations constructed from adaptive high-regularity truncated hierarchical splines, as employed in the IGA paradigm. 27,28 Our choice for high-regularity splines is motivated by the superior approximation properties of such spline functions on a per-degree-of-freedom basis 29 relative to traditional piecewise-polynomial finite-element approximations. The considered adaptive strategy comprises a two-level hierarchical a posteriori error estimate, that is, the error associated with each mesh in the sequence of adaptive refinements is estimated by comparing the approximation to the approximation on a mesh in which all elements are uniformly subdivided. An important advantage of this approach is that it is semi-automatic, in the sense that it does not require the derivation of an error-residual relation. In contrast, a conventional residual-based error estimate would insist on the derivation of a sharp upper bound on the element-wise approximation error in terms of the local residuals, which is prohibitively complicated in elasto-capillary FSI on account of the multitudinous fields that appear in the formulation and the complex form of the equations, involving various intricate interface and coupling terms. The hierarchical a-posteriori error estimate directs a support-based refinement procedure. 30 To avoid the nonrobustness of partitioned iterative solution methods for elasto-capillary FSI, 21 we apply a monolithic solution procedure to solve the aggregated binary-fluid-solid-interaction (BFSI) problems in each time step and each iteration of the adaptive-refinement process. To enhance the robustness of the monolithic solution procedure, we introduce an -continuation process, that is, the value of the interface-thickness parameter, , is artificially enlarged on the coarsest meshes in the adaptive procedure.
To assess the properties of the presented adaptive IGA methodology for elasto-capillary FSI, we conduct numerical computations for a configuration pertaining to a droplet on a soft solid substrate, and we present a comparison to experimental data from Reference 31. The remainder of this article is structured as follows. Section 2 presents the BFSI problem. In Section 2.4 we elaborate on the error-estimation method and the adaptive-refinement process. This section also presents details of the applied locally refined truncated hierarchical B-spline approximations and of the solution procedure. Section 4 presents numerical experiments and results. Section 5 is concerned with concluding remarks.

ELASTO-CAPILLARY FSI MODEL
The elasto-capillary fluid-solid-interaction model considered in this work is similar to that in Reference 21 with the exception of the NSCH model, which is fully incompressible in accordance with Reference 15 as opposed to the quasi-incompressible model in Reference 21. The main aspects of the model are repeated here for completeness and coherence of the presentation. We refer to Reference 21 for a detailed derivation and discussion of the model; see also Reference 26 for an elasto-capillary FSI model based on the Navier-Stokes-Korteweg equations. We consider a binary fluid composed of two immiscible constituents separated by a thin-but-finite interface; see Figure 1. The binary fluid is contiguous with a deformable elastic solid. The BFSI problem is set on a time interval (0, t fin ) ⊆ R >0 and two simply-connected time-dependent open subsets Ω f t ⊂ R d (d = 2, 3) and Ω s t ⊂ R d , accommodating the binary fluid and solid, respectively. The solid-fluid interface Γ t ∶= Ω f t ∩ Ω s t corresponds to the intersection between the boundaries of the fluid and solid domains. We assume that the time-dependent configuration Ω t ∶= int(cl Ω f t ∪ cl Ω s t ) (where int and cl denote the interior and closure, respectively) can be viewed as the image of a continuous transformation T t acting on the initial configuration Figure 1. Note that, in particular, the transformation T t is continuous across the solid-fluid interface. Moreover, we assume that the restriction of T t to both Ω f 0 and Ω s 0 is a diffeomorphism.

Diffuse-Interface Binary-Fluid Model
The components of the binary fluid are equipresent according to the volume fractions l , a ∶ (0, t fin ) × Ω f t → [0, 1]. The volume fractions form a partition of unity. The presence of each constituent, locally in space and time, is identified by the order parameter: and l = (1 + )∕2 and a = (1 − )∕2. The pure l-species and a-species phases are indicated by = 1 and = −1, respectively. Intermediate values of the order parameter indicate a mixture. Denoting by l and a the specific densities of the two constituents, the mixture density ∶ (0, t fin ) × Ω f t → R >0 is related to the volume fractions by = l l + a a and hence: We assume that l ≥ a > 0. Denoting by u ∶ (0, t fin ) × Ω f t → R d the volume-averaged mixture velocity, we consider the incompressible NSCH diffuse-interface model: 15 with a relative mass flux 15 and Ψ( ) = 1 4 ( 2 − 1) 2 the standard double-well potential for the energy density of the mixture. (3) represent fluid pressure and chemical potential, respectively. We assume that the mobility, m > 0, is constant.
The fluid stress in (3) comprises a viscous component, , a capillary component, , and the isotropic component pI. The constitutive relations for the viscous stress and capillary stress are: The mixture viscosity ∶= ( ) is defined via the Arrhenius relation: with Λ = l m a ∕ a m l as the intrinsic liquid-ambient volume ratio, that is, the ratio of the molar volume of the liquid to the molar volume of the ambient fluid, and m a and m l the molar densities. The parameter is related to the liquid-ambient surface tension la by 2 √ 2 = 3 la . The parameter > 0 controls the thickness of the diffuse interface; see for instance Reference 21. The trace of the capillary stress tensor in (6) is generally non zero and, hence, the pressure p does not admit the usual interpretation as the isotropic part of the stress tensor.
Remark 1. The order-parameter, , admits an unambiguous physical interpretation according to (1), provided that it only assumes values in the interval [−1, 1], almost everywhere. However, in general, this condition on the range of must be dismissed to ensure existence of solutions to the NSCH system; see, for instance, Reference 32. This implies that can exceed the bounds ±1. Accordingly, the interpolated density (2) may then attain nonpositive values, in particular if the density contrast is large. Motivated by the results in Reference 32, to avoid such degeneration of the interpolated mixture density we extend the mass-density interpolation into | | ≥ 1 by with regularization parameter = a ∕( l − a ). The mixture-density function (8) coincides with the C 1 -continuous soft-clipped linear interpolation introduced in Reference 33.
Remark 2. The mixture viscosity in NSCH models is generally also interpolated by means of a (clipped) linear function analogous to (8). However, whereas the linear interpolation of the mixture-mass density in (2) is imposed by conservation of mass, the interpolation of viscosity is merely a constitutive model, and alternate models can be considered. Unlike the linear interpolation of viscosity, the Arrhenius mixture-viscosity model (7) considered in this work has a rational basis in physical chemistry via the Grunberg-Nissan mixing rule. Moreover, if the the interpolation is based on the volume fraction (Λ = 1), the mixture-viscosity relation (7) is nondegenerate, and the mixture viscosity is bounded away from zero, that is, ( ) ≥ c > 0 for some constant c. In this work we will restrict ourselves to the volume-fraction-based Arrhenius relation.
A comprehensive treatment of initial and boundary conditions for (3) is beyond the scope of this work; see References 21,22 for some further considerations. The boundary conditions at the solid-fluid interface are addressed in Section 2.3.

Incompressible hyperelastic solid model
The deformation of the soft solid substrate is described by an incompressible hyperelastic-solid model. We use the( ⋅) symbol to distinguish functions that are defined in the reference configuration. Denoting the mass density of the solid byr ∶ Ω s 0 → R >0 , the solid deformation̂∶ (0, t fin ) × Ω s 0 → Ω s t satisfies the equation of motion and the incompressibility constraint: withF = Grad̂the deformation tensor andP =Q −̂I the first Piola-Kirchhoff stress tensor. Furthermore, Div(⋅) and Grad(⋅) denote the divergence and gradient operators in the reference configuration, respectively, and̂∶ (0, t fin ) × Ω s 0 → R represents solid pressure. The constitutive relation for the first Piola-Kirchhoff stress follows from the derivative of the strain-energy density to the deformation tensor.
Because the reference configuration coincides with the initial configuration, the initial condition̂(0, ⋅) = (⋅) holds. Boundary conditions at the solid-fluid interface are presented in Section 2.3. Auxiliary initial and boundary conditions for (9) are not further considered here.

Interface conditions
Similar to conventional FSI problems, the fluid and solid subproblems in BFSI problems are interconnected via a topological condition, a kinematic condition and a dynamic condition. In addition, the binary-fluid satisfies a wetting condition, which expresses the affinity of the solid substrate to the two fluid components, and a boundary condition imposing the impermeability of the solid to the fluid species. The topological and kinematic conditions are identical to those in conventional FSI. The dynamic condition deviates from the conventional case, owing to the contribution of the capillary-stress tensor, , to the fluid traction and the surface tension carried by the solid-fluid interface, which leads to a traction discontinuity across the interface in accordance with the Young-Laplace relation.
The topological condition implies that the fluid and solid subdomains are contiguous at the interface. This condition is trivially satisfied by the continuity of the transformation map T t ∶ Ω 0 → Ω t . We construct the transformation map as follows: where Êis an extension of the solid deformation,̂, onto the fluid domain. The extension operator corresponds to the deformation of a linear-elastic pseudo-solid. 34 The kinematic condition imposes continuity of the binary-fluid and solid velocities at the interface: where T # t denotes the pullback (by T t ) to the reference configuration, that is, (T # t u)(X) = u(T t X). The dynamic condition imposes equilibrium of the fluid and solid tractions and the traction exerted on the interface by solid-fluid surface tension. Solid-fluid surface tension can generally be ignored in conventional FSI problems, but is of fundamental importance in elasto-capillary FSI for soft solids. 21 The traction induced by the solid-fluid surface tension satisfies the Young-Laplace relation for nonuniform surface tension: with sf ∶= sf ( ) the composition-dependent solid-fluid surface tension (see (16) below), the additive curvature of the interface, and grad Γ (⋅) the tangential gradient on Γ t ; see, for example, References 21,22,35. We adhere to the convention that curvature is negative if the center of the osculating circle in the normal plane is located in the fluid domain. The binary fluid in (3) exerts a traction on the interface according to: The binary-fluid traction contains a capillary contribution in addition to the usual pressure and viscous traction. The dynamic interface condition imposes that: with t s the traction exerted by the solid substrate on the interface. By Nanson's formula it holds that T # t t s = −Pn s 0 , where corresponds to the ratio of surface measures in the actual and reference configuration.
The wetting condition represents the affinity of the two fluid components to the solid substrate at Γ t : where sf represents solid-fluid surface tension according to with sl ≥ 0 and sa ≥ 0 the solid-fluid surface tensions of the liquid species and ambient species, respectively; see also References 21,22,36-38. Condition (15) imposes that in equilibrium the diffuse fluid-fluid meniscus intersects the solid-fluid interface at the equilibrium contact angle e = arccos(( sa − sl )∕ la ); see also References 21,38. It is to be mentioned that the equilibrium contact angle is defined in the usual manner as the angle (interior to the liquid) between the meniscus of a sessile liquid droplet and a rigid, planar solid surface that supports it. The impermeability condition insists that the solid surface is closed to the two fluid components. Noting that m n corresponds to the diffusive volume-fraction flux, the impermeability condition implies n = 0.

Weak Formulation of the BFSI Problem
This section presents a weak formulation of the BFSI problem. A comprehensive derivation of the weak formulation is not envisaged here, and we refer to Reference 21 for further details. In the weak formulation, we consider a global deformation field̂∶ (0, t fin ) × Ω → R d . Its restriction to Ω s 0 represents the solid deformation, while its restriction to Ω f 0 represents the deformation of the fluid domain. Similarly, we define a global velocity field,v ∶ (0, t fin ) × Ω → R d , which coincides with the domain deformation velocity,v = t̂. The fluid velocity, u, is expressed in the reference domain and with respect to the background velocity,v, as T # t u =v +û. The kinematic condition (11) is fulfilled by insisting that u| Γ 0 = 0. The dynamic condition is incorporated in the weak formulation by applying a global testfunctionŷ ∶ Ω → R d for the equations of motion of the fluid and the solid 21 in (3 1 ) and (9 1 ), respectively. The motion of the fluid domain is accounted for by a weak ALE formulation of the NSCH system (3).
To provide a setting for the weak formulation, we define aggregated ambient trial and test spaces. For conciseness, we will just introduce the basic structure of the ambient spaces; see Reference 21 for further discussion. The test and trial functions are defined on the reference domain and are transported to the actual (time-dependent) domain in the weak formulation as appropriate. We set The pullback of a function from the reference configuration to the actual configuration by the inverse transformation T −1 t will be indicated by omitting the( ⋅) symbol, that is, . The test space for the weak formulation of the BFSI problem isÛ. We note that the trial and test spaces are generally subject to auxiliary conditions, for example, related to essential boundary conditions, but we will not consider such conditions further.
The weak formulation of the aggregated BFSI problem can be condensed into: where l ∶Û → R denotes a linear form containing exogenous data and the semilinear forms a 0 , a 1 ∶Û ×Û → R are defined by with (⋅, ⋅) a suitable inner product. The tensor̂encodes the constitutive relation for the pseudo-solid problem that provides the deformation of the fluid domain. The semi-linear form a 0 represents the time-derivative terms. The first three terms in (19b) represent the equation of motion for the solid, including the weak formulation of the load induced by solid-fluid surface tension. 21 The fourth term imposes the pseudo-solid behavior on the deformation in the fluid domain, | Ω f 0 . Terms 5 to 8 represent the ALE formulation of the equation of motion of the binary fluid in (3 1 ). The ninth term represents the ALE formulation of the spatial part of the convective Cahn-Hilliard equations in (3 3 ). Terms 10 to 12 correspond to the definition of the chemical potential in (3 4 ). Term 13 results from substituting n = − ′ sf ( ) and serves to impose the wetting condition (15). The penultimate term in (19b) enforces solenoidality of the mixture velocity in accordance with (3 2 ). The ultimate term in (19b), in conjunction with the ultimate term in (19a), imposesv = t̂.

ERROR ESTIMATION AND ADAPTIVE APPROXIMATION
We consider adaptive approximations of (18) based on Truncated Hierarchical B-Splines (THB-Splines 39 ). The local refinements are directed by a hierarchical error estimate following the standard recursive SEMR (Solve → Estimate → Mark → Refine) process. 40,41 Motivated by the nonrobustness of partitioned iterative solution methods for elasto-capillary FSI 21 problems, we apply a monolithic solution procedure with continuation to solve the aggregated BFSI problems in each time step and each iteration of the adaptive-refinement process.

Truncated hierarchical B-Splines
THB-splines can be constructed from unions and linear combinations of hierarchical B-splines. Multivariate B-splines are formed as tensor products of univariate B-splines. To define the elementary univariate B-spline basis, we consider a partition of the unit interval (0, 1) by a sequence of nondecreasing coordinates (or knots), Ξ = { 1 , 2 , … , k+m+1 } with 1 = 0 and k+m+1 = 1, where k and m denotes the order and dimension of the B-spline basis, respectively. The knot vector Ξ serves as a substructure for the B-spline basis, comprised of the piece-wise polynomials N i,k , recursively defined via the Cox-De Boor relation: 42,43 For j ≥ 1, the B-splines are supported on multiple knot spans (cf. elements), that is, the line segment between consecutive distinct knots in Ξ, and the support increases with order j. If all knots in Ξ have multiplicity one, that is, without knot repetition, the B-splines according to (20) are C k−1 continuous. If a knot i is repeated with multiplicity a i , the regularity of the corresponding B-spline reduces locally to C k−a i . Generally, the first and final knots, 1 and k+m+1 , are assigned multiplicity k + 1, rendering the basis interpolatory at the boundaries, {0, 1}, that is, 1 A refinement of a B-spline basis can be obtained by the process of knot insertion. Considering a knot vector Ξ 0 and the augmented vector Ξ 1 obtained by adding the midpoints of distinct knots in Ξ 0 , that is, by bisecting the knot spans of Ξ 0 , it holds that S(Ξ 1 ) ⊃ S(Ξ 0 ). This process of knot insertion can be repeated to generate a sequence of knot vectors Ξ 0 , Ξ 1 , Ξ 2 , … and a corresponding sequence of nested B-spline spaces S(Ξ 0 ) ⊂ S(Ξ 1 ) ⊂ S(Ξ 2 ) ⊂ …. The knot-insertion construction extends directly to tensor products of knot vectors and can hence serve to define a hierarchy of nested multivariate B-spline spaces.
To facilitate the presentation, in the definition of locally refined hierarchical d-variate B-spline spaces, we restrict ourselves here to isotropic and uniform spline spaces on the unit cube, (0, 1) d . Given a knot-span length h 0 ∈ 1∕N, spline order k ∈ N and regularity index ∈ Z ≥0 , we denote by Ξ ( ∈ Z ≥0 ) the family of open knot vectors: by B k, the univariate B-spline basis B(Ξ ) and by S k, (0, 1) = S(Ξ ) the corresponding space. The affixes k, will be tacitly suppressed if these are clear from the context or irrelevant for the exposition. Let us note that by virtue of the multiplicity k − of the interior knots, S k, (0, 1) ⊂ C (0, 1). The d-variate B-spline basis at level corresponds to the tensor product B d , defined as the collection of d-products of the univariate basis functions. The d-dimensional level-knot spans are defined as the collection Denoting by L the maximum refinement level, we consider a nested sequence of nonempty subdomains (refinement regions), D 0 , … , D L , and its completion by empty sets: subject to the condition that each subdomain D corresponds to the union of a collection of level-knot spans: with 2 the powerset of . Given the uniform spline bases B d and the refinement regions D , we associate to each level a local refinement L : where supp • ( ) denotes the open support of . In particular, L comprises all level-B-splines whose supports are contained in the level-refinement region and are not completely covered by the level-( + 1) refinement region. The locally refined HB-basis is then defined as the union ∪ L =0 L . The truncation of Hierachical B-splines is motivated by the fact that Hierarchical B-spline bases do not generally satisfy the partition-of-unity property and may exhibit excessively many overlaps between basis-function supports. 39 The truncation operation is based on the observation that by virtue of the nesting property of spline spaces, each level-B-spline can be expressed as a linear combination of level-( + 1) B-splines. Hence, each B-spline ∈ L L−1 can be expressed as a linear combination: The truncation consists in suppression of the part of N L−1 that is already contained in L L , that is, by setting the coefficient c 0 in (26) to zero: Note that if supp • ( ) ∩ supp • ( ) = ∅ for all ∈ L L , then tr = . The level L − 1 truncated local refinement can then be formed as L tr L−1 = { tr ∶ ∈ L L−1 } and it holds that (spanL L ⊕ spanL L−1 ) = (spanL L ⊕ spanL tr L−1 ). Similarly, each spline ∈ L L−2 can be expressed as a linear combination: for certain c 0 , c 1 , generally distinct from c 0 , c 1 in (26). The part contained in L L ∪ L tr L−1 can again be suppressed, leading to tr . Applying this procedure for all ∈ L L−2 yields the truncated local refinements L tr L−2 . Continuing this representation and truncation process recursively from level L to level 0 yields a truncated local refinement at each level. The THB-spline basis is then obtained as the union of the local refinements over all levels: Remark 3. In Reference 39, the hierarchical truncation procedure is executed in a coarse-to-fine arrangement. We opt for the aforementioned fine-to-coarse arrangement, progressing from level L to level 0. Because the THB bases in both arrangements span the same space and both form a partition of unity, the bases are in fact identical. There can, however, be differences in the representation of each basis function in the THB basis with respect to the underlying B-splines.

Approximation of the BFSI Problem
We consider discrete approximations of (18) based on an implicit Euler finite-difference approximation in the temporal dependence and IGA approximation in the spatial dependence based on the THB-spline spaces introduced in Section 3.1.
We reintroduce the affixes k, of the hierarchical space H k, L to indicate the order and regularity of the considered spline spaces. Let T f ∶ Ω f 0 → (0, 1) d and T s ∶ Ω s 0 → (0, 1) d respectively represent transformations from the fluid and solid reference domains to the unit cube. The IGA approximation space for the BFSI problem comprises (cf.Û in (17)): at the solid-fluid interface coincide. This condition is implicit in the construction of the approximation space for the deformation,̂, and the corresponding velocity,v, which is composed of the THB-spline spaces in the fluid and the solid, and is supposed to be continuous across the interface.
Remark 4. For k ≥ 3, the deformation (and corresponding velocity) in the fluid and solid domains are more regular than C 0 -continuous, and the reduction to C 0 continuity at the solid-fluid interface in (30) requires repetition of the knots at the interface. This reduction of the regularity at the interface is consistent with the fact that the deformation in the solid and fluid subsystems is generally described by distinct models, that is, an incompressible nonlinear hyper-elastic model in the solid subsystem and a compressible linearly elastic pseudo-solid model in the fluid subsystem, and hence the deformation is generally only continuous across the interface, but not more regular. Numerical experiments indicate that retaining C k−2 -continuity of the deformation across the interface generally causes non-robustness in the Newton iteration (results not displayed).
Remark 5. The fluid velocity and solid deformation are approximated with B-splines of order k ≥ 2 and regularity k − 2, while the pressure in the fluid,p L , and in the solid,̂L, are approximated with B-splines of order k − 1 and regularity k − 2. These velocity/pressure and deformation/pressure pairs correspond to IGA extensions of the Taylor-Hood family of finite-element pairs. 44 For k = 2, the standard C 0 -continuous Q 2 ∕Q 1 Taylor-Hood finite-element pair is recovered.
Exploiting the rotational symmetry, we furnish Ω with an approximation space that conforms to (30) in the radial-longitudinal dependence and that is constant in the azimuthal dependence. Denoting byÛ L a two-dimensional approximation space conforming to (30), the rotationally symmetric approximation space on the cylinder,Û cyl L , corresponds to the pullback by T −1 cyl of the constantly extended spaceÛ L : In the presentation of the approximation of the BFSI problem (18), we will restrict ourselves to the general case, but the formulation extends directly to the rotationally symmetric case by replacingÛ L byÛ cyl L . To discretize (18) in time, we introduce a uniform partition of the time interval under consideration, (0, t fin ), into subintervals of length . We denote byÛ n L ∈Û n L an approximation to the solution,Û(t n ), of (18) at time t n = n . It is to be noted that the approximation space,Û n L , generally depends on the index n. The approximation spaces {Û n L } are not predefined, but are constructed in each time step via an adaptive-refinement process; see Section 3.3. The discrete formulation is based on a backward-Euler approximation in the time dependence. In addition, the semilinear form in the discrete formulation is consistently modified with respect to a 1 in (19b) by a decomposition of the double-well potential 18,45,46 and a skew-symmetric formulation of the convective form, 47 to enhance the stability of the formulation. Given a suitable initial elementÛ 0 L ∈Û 0 L , for instance obtained by a projection of the initial data, the approximation of the BFSI problem can be condensed into: FindÛ n L ∈Û n L such that: for n = 1, 2, …. In the discrete approximation (33), we take the inner product that appears in the semi-linear forms a 0 , a 1 in (19) as the L 2 (Ω, R d )-inner product, that is, (v,ẑ) = ∫ Ωv ⋅ẑ. The semi-linear form s 0 in (33) includes the expansive part Ψ + ( ) = 1 4 4 − 3 2 2 of the double-well potential: Accordingly, in (33) the contractive part of the double-well potential, Ψ − = Ψ − Ψ + , is treated implicitly, while the expansive part Ψ + is treated explicitly. This decomposition of the double-well potential serves to enhance the stability of the formulation. 18,45,46 The semi-linear form s 1 is defined as: One can infer that s 1 constitutes a partition of zero if v + u is solenoidal. Hence, s 1 is consistency presevering. Noting that the first term in s 1 coincides with the convective form in (19b), addition of s 1 in (33) renders the convective part of the formulation skew-symmetric, in the sense that the volumetric terms, except for the term involving grad , cancel if y is replaced by v + u. This modification of the semi-linear form serves to cancel artificial energy production via the equation-of-motion due to inexact solenoidality in the Taylor-Hood approximation. 47

Error Estimation and Adaptive Refinement
We consider an adaptive-refinement procedure directed by an a posteriori error estimate. Motivated by the complexity of the BFSI problem, notably the fact that it involves seven different fields, we revert to a generic two-level hierarchical error estimate, in which the difference between the approximation inÛ L and its so-called superrefinement is used to estimate the error for each field. An upper bound of the error estimate is subsequently decomposed into support-wise error contributions, which are used to mark the function supports that need to be refined to most effectively reduce the error. To elaborate the two-level hierarchical error-estimation procedure, we consider (33) in a time step n ∈ {1, 2, …}, and we assume that the approximations in all previous time steps 0, … , n − 1 have been satisfactorily resolved. We conceive of the approximation spaceÛ n L as the currently highest level of refinement in the adaptive procedure in time step n, and we denote byÛ n L the corresponding approximation. Denoting byÛ n L the superrefinement ofÛ n L , that is, the approximation space that is obtained by uniform subdivision of each knot span ofÛ n L into 2 d equal parts (see Figure 2), and byÛ n L the corresponding approximation of (33), the error inÛ n L can be estimated asÛ n L −Û n L . To use the error estimateÛ n L −Û n L to direct an adaptive-refinement process, the error estimate must be decomposed into local contributions that can serve to introduce the refinements. In adaptive-refinement procedures for conventional finite-element methods, the local error contributions are generally assigned to elements. In IGA approximations, it is more convenient to assign the error contributions to basis-function supports. 30,48,49 To measure the error in the various field ofÛ n L , we introduce the norm: F I G U R E 2 Illustration of the adaptive refinement procedure for three fields. The red line in each panel represents the solid-fluid interface, but the interface does not appear explicitly in the refinement procedure. Panel (A) depicts the original refinement regionŝ D 0 = Ω (light grey) andD 1 (medium grey) and knot spans corresponding to the hierarchical basis. Panel (B) displays the knot spans of the super-refinement, and the function supports marked for refinement by field 1 (blue), field 2 (green) and field 3 (yellow). Multicolored regions have been marked by multiple fields. Panel (C) indicates the extended refinement regionD * 1 (medium grey) and the new level 2 refinement regionD * 2 (dark grey). The crosses indicate auxiliary refinement that are introduced by the completion procedure. Panel (D) displays the knot spans of the refined hierarchical basis, and the extended and completed refinement regionsD 0 (light grey),D 1 (medium grey) andD 2 (dark grey).
To each field f ∈ {û,p,̂,̂,̂,v,̂} =∶ F and each corresponding basis function in the super-refinement, we assign a support-wise error indicator according to: where || ⋅ || supp • ( ) denotes the restriction of the norm pertinent to the field f to the open support of , and A f L denotes the basis of the super-refinement for the field f . For instance, By virtue of the triangle inequality and the overlap of the basis-function supports, the support-wise error indicators provide an upper bound to the error estimate: The objective of the adaptive-refinement procedure is to effectively reduce the error by refining basis functions associated with the principal error contributions. In particular, for each field f ∈ F (or a suitable subselection of fields; see Remark 7) we mark the largest error contributions according to a Dörfler-type marking, 41 that is, we select a minimal set of basis for some envisaged contraction coefficient a f ∈ (0, 1). It is to be noted that according to (40), if the adaptive refinement indeed completely eliminates the error contributions associated with ∈ M f , then the (upper bound to the) error in the considered field reduces by a factor a f . The refinements are introduced by extending the refinement regions that form the support structure for the locally refined hierarchical basis; cf. Section 3.1 and see Figure 2 for an illustration. We denote by the refinement regions corresponding to the hierarchical basisÛ n L in the reference configuration. The extension of the refinement regions can be conveniently separated into two operations. Noting that by virtue of the construction of THB bases, each basis function ∈ M f has an intrinsic hierarchical rank, rank( ) ∈ {1, 2, … , L + 1}, we first add the support of marked basis functions with hierarchical rank to the refinement domainD : where supp c ( ) denotes the closed support of . Second, we perform a completion step, by incorporating all level-( − 1) knot spans that intersect withD * , to formD : where −1 denotes the set of level-( − 1) knot spans in the reference configuration. The unary operator ← indicates an update of the left member. The updated refinement domains {D } L+1 =1 are subsequently used to construct the refined THB-spline basis, following the process described in Section 3.1. The refinement process is shown in Figure 2.
Remark 6. In principle, the refinement regionsD * could be used directly as level refinement domains. The extension ofD * to level-( − 1) knot spans in (43) serves to retain a regular structure of the knot spans in the THB-spline basis to carry out numerical quadrature. In particular, the collection of all knot spans for the locally refined THB-spline basis, provides a shape-regular tessellation of the domain Ω.
Remark 7. In the exposition of the adaptive-refinement process, we have assumed that all BFSI-fields are included in the error estimate and in the marking-and-refinement process. However, in some cases, it is appropriate to restrict the refinement process to a subselection of the fields or, equivalently, to set f = 0 for certain fields. In particular, if one is interested in equilibrium configurations 21 of the BFSI problem, then the velocity fields,v andû, and the chemical potential,̂, generally vanish. For equilibrium and near-equilibrium configurations, these fields can therefore be excluded from the marking procedure.

Monolithic solution procedure and -continuation
In each time step and each iteration of the adaptive process, two nonlinear discretized BFSI problems must be solved, namely (33) for the hierarchically refined space,Û n L , and for the corresponding superrefinement,Û n L . Solution methods for algebraic systems deriving from FSI problems can generally be classified into partitioned methods and monolithic methods. Monolithic methods solve the aggregated FSI problem in a fully coupled manner. Partitioned methods generally rely on a decomposition of the algebraic system into fluid and solid subsystems, which are then solved alternatingly until a prescribed convergence criterium is satisfied. The convergence behavior of the basic partitioned method for FSI, namely subiteration with under-relaxation, has been investigated for elasto-capillary FSI in Reference 21. The results in Reference 21 indicate severe nonrobustness of partitioned methods for elasto-capillary FSI.
Motivated by the nonrobustness of partitioned methods for elasto-capillary problems, in this work we apply a monolithic solution procedure. The monolithic method consists of a Newton method enhanced with line search, in combination with direct solution of the linear tangent problems in the Newton method. The tangent problems are based on a consistent linearization of the semi-linear form −1 a 0 + a 1 + s 0 + s 1 in (33). The required Fréchet derivatives of the complete semi-linear forms in (19) are determined via automatic differentiation. It is to be noted that manual differentiation of (19) to, in particular, the deformation is profoundly complex, due to the implicit dependence of Ω f t and the surface gradients in the solid-fluid-surface-tension term on the deformation. Automatic differentiation extracts the Fréchet derivative directly from the implementation of the semi-linear form itself.
Without further precautions, in the initial stages of the adaptive-refinement process the minimal mesh width (knot span) is generally significantly larger than the thickness of the diffuse interface. Such severe under-resolution of the diffuse interface can impede convergence of the Newton procedure. To enhance the robustness of the monolithic solution procedure, we apply an -continuation process, that is, on the coarsest levels in the adaptive procedure the value of the interface-thickness parameter, , is artificially enlarged. Specifically, we select an initial interface-thickness parameter 0 = 2 K such that 0 is close to the baseline mesh width h 0 and K is sufficiently small relative to the maximum number of refinement steps, L max . At each level L = 0, 1, 2, … , L max of the iterative refinement process, instead of solving (33) with the targeted interface-thickness parameter , we solve (33) with an interface-thickness parameter L = max(2 K−L , 1). Hence, in the initial K − 1 stages of the adaptive-refinement procedure, the interface thickness is halved in each iteration, while in the final L max + 1 − K stages the interface thickness is fixed at its prescribed value .
Remark 8. The -continuation process amounts to a modification of the adaptive-refinement process described in Section 3.3, in that the error estimate that directs the adaptive procedure is based on L instead of in the first K − 1 refinement stages. This modification relies on the assumption that the refinement strategy is insensitive to the precise value of the interface-thickness parameter if the diffuse interface is not adequately resolved. However, this assumption appears not more restrictive than the intrinsic assumption that the a posteriori error estimate and the derived error indicators can effectively direct the adaptive-refinement process even if the approximation is severely under-resolved.

NUMERICAL EXPERIMENTS
To assess the properties of the adaptive IGA framework for elasto-capillary FSI outlined in Sections 2 and 3, we conduct numerical experiments pertaining to a sessile droplet on a soft solid substrate; cf, Figure 1. We provide a comparison of equilibrium solutions of the BFSI problem obtained with the adaptive IGA simulation framework with experimental data. 31 In addition, we consider transient results and aspects of the adaptive-refinement process.

Equilibrium solutions
We regard a liquid droplet set in an ambient fluid on a soft solid substrate, see Figure 3. The corresponding actual experiment as conducted in Reference 31 pertains to a glycerol droplet in air on a silicone-gel layer. The numerical setup is rotationally symmetric. The initial configuration of the binary-fluid domain corresponds to a circular cylinder with radius 500 m and height 500 m. The initial configuration of the solid domain is a circular cylinder with radius 500 m and height 50 m. The boundary conditions of the BFSI problem at the external boundaries are shown in Figure 3. The solid substrate is secured at the bottom boundary and guided at the lateral boundary. The pseudo solid is guided at its external boundaries. The boundary conditions at the lateral and longitudinal boundary for the binary fluid are shown in Figure 3. The initial condition for the order parameter corresponds to a spherical-cap-shaped liquid volume, positioned such that the fluid-fluid meniscus intersects the undeformed solid substrate at the equilibrium contact angle, in particular, with R l the radius of the spherical cap, x c = −(R l cos e ) e 3 its center, and e 3 the unit vector in longitudinal direction. 26 The droplet volume, 13.7 nL, has been selected such that the radius of the circular footprint corresponds to the value of 176.7 m reported in Reference 31. The surface tension of the fluid-fluid meniscus is la = 46 mJ∕m 2 . The solid-fluid surface tension of the liquid (resp. ambient fluid) is sl = 36 mJ∕m 2 (resp. sa = 31 mJ∕m 2 ). The shear modulus of the solid substrate is set to G = 1 kPa in accordance with the specification of the Young's modulus in Reference 31 and the incompressibility of the material. We set the interface length-scale parameter to = 1 m, which is approximately 200× smaller than the droplet radius and approximately 50× smaller than the elasto-capillary length, la ∕G. The substrate is modeled as an incompressible hyperelastic solid with a neo-Hookean strain-energy-density function: withÊ ∶=Ê(F) = 1 2 (F TF − I) the Green-Lagrange strain tensor. The neo-Hookean constitutive relation (46) provides a more appropriate model for silicon gel than the St. Venant-Kirchhoff relation considered in Reference 21. The pseudo-solid is equipped with linearly elastic constitutive behavior, with unit modulus of elasticity and Poisson ratio 1∕4.
Our primary interest pertains to equilibrium solutions of the elasto-capillary FSI problem. The material mass densities, the viscosities of the liquid and ambient fluid, and the mobility only affect the dynamic behavior, but these are inconsequential for the equilibrium configuration. The liquid (glycerol) and ambient fluid (air) considered in the actual experiment have density and viscosity 1.26 kg∕L and 1.41 Pas, and 10 −3 kg∕L 3 and 20 Pas, respectively. The density and viscosity ratio of the liquid and ambient fluid exceed 10 3 and 10 5 , respectively, which compromises the stability of the numerical approximation for small values of and large values of . To enhance the stability of the computation, the density and viscosity of the ambient fluid are matched to those of the liquid. An additional advantage of this matched-density, matched-viscosity setting is that the obtained results are essentially independent of the applied NSCH model, as all NSCH models reduce to the conventional model-H 14 in this case. For the solid, we select a mass density of 1 kg∕L, appropriate to silicone gel. A synopsis of the parameter settings for the considered test case is shown in Table 1. We consider numerical results for k = 3 splines; cf, (30). In particular, the deformation and domain velocity in the fluid and solid domains are approximated with C 1 -continuous cubic splines, the relative fluid velocity is also approximated with C 1 -continuous cubic splines, the pressure fields in the fluid and the solid are approximate with C 1 -continuous quadratics, and the order parameter and chemical potential are approximated with C 2 -continuous cubic splines. The time-step size is set to = 10 3 s. In each time step, the fluid domain is initially covered with a uniform mesh comprising 20 × 20 elements in the radial/longitudinal dependence, and the solid domain is covered with a uniform mesh of 20 × 2 elements. The mesh width of the initial mesh is therefore h 0 = 25 m. In each time step, L max = 10 refinement steps are executed. Accordingly, the smallest elements in the adaptively refined mesh are up to 2 10 smaller than the initial mesh width. In this regard, it is important to note that the maximum hierarchical level in the refinement can be smaller than the number of refinement iterations, because the set of basis functions selected for refinement, ∪ f ∈F M f , does not necessarily contain basis functions at the highest hierarchical level. Considering that our interest is essentially restricted to equilibrium configurations, the marking procedure only includes the order parameter, the fluid and solid pressure fields, and the deformation; see Remark 7. The coefficients in the Dörfler-type marking are set tô= 0.5,̂p = 0.9,̂= 0.5 and̂= 0.5. The relatively large value of the refinement factor̂p is motivated by the fact that we anticipate that the error localizes on the diffuse interface and that, in principle, all elements on the interface need to be refined to achieve an effective error reduction. -continuation is applied in the first K = 3 stages of the refinement procedure. Figure 4 presents a comparison of the computed near-equilibrium solid-fluid surface deformation at 15 ms with experimental data. The computed result exhibits excellent agreement with the experimental data, especially in the wetting-ridge region. The wetting-ridge height closely matches the experimental data. The indentation below the droplet, caused by the elevated pressure in the droplet and incompressibility of the solid, is slightly deeper than for the experimental data. The present results exhibit much better agreement with the experimental data than those presented in References 21,22,26. This improvement in the results is predominantly caused by the modified soft-solid constitutive behavior, namely, neo-Hookean in the present work vs St. Venant-Kirchhoff in the aforementioned references. Figure 5 presents an overview of the near-equilibrium solution of the BFSI problem at t = 15ms, including a magnification of the contact-line region. In the overview, the mesh in the binary-fluid domain has been suppressed to provide a clear view of the diffuse fluid-fluid interface. The error-estimate-based adaptive-refinement procedure automatically concentrates the refinements in the vicinity of the diffuse interface and near the apex of the so-called wetting ridge. By virtue of the higher-order approximation provided by the C 1 splines, the diffuse interface is essentially resolved with 6 levels of refinement, corresponding to a mesh width of approximately h 6 = 0.39 m ≈ ∕2. Indeed, after L max = 10 refinement steps, the diffuse interface is covered for the largest part with splines of hierarchical level 6, and the maximum hierarchical level 7 only occurs in isolated locations. The colors in five indicate the pressure, that is, the scaled trace of the binary-fluid stress tensor, p − tr ∕d. The pressure in the binary fluid exhibits a strongly localized minimum at the diffuse fluid-fluid interface. In the sharp-interface limit → +0, the traction exerted by the binary fluid F I G U R E 5 Near-equilibrium binary-fluid-solid interaction solution for test case 1 at t = 15ms. The inset displays a magnification (25×) of the contact-line region, including the adaptive mesh in the binary-fluid. Colors represent the pressure at the contact line reduces to a delta distribution, and the solid-fluid interface exhibits a sharp kink at the apex. 21,26 For the considered value of the interface-thickness parameter, = 1 m, the interface deformation at the apex exhibits a small but finite radius of curvature. From Figure 5 one can observe that the pressure in the droplet is virtually uniform at an elevated level of approximately 520 Pa relative to the ambient pressure. This elevated pressure is close to the theoretical Laplace pressure 2 la ∕R in a droplet on a rigid substrate with radius R = 178 m and surface tension la = 46 mN∕m.

Evolution of the Wetting Ridge
We next consider the dynamics of the aggregated elasto-capillary FSI problem. We regard a configuration with parameter values according to test case 2 in Table 1. The settings for the solid substrate, the liquid and the ambient-fluid are representative of the materials considered in the experimental setup in Reference 31, except for the ambient-fluid viscosity, which has been increased to 0.1 Pas to retain stability in the numerical computation. To monitor the dynamics of the wetting ridge, the time step has been set to = 10 s. Figure 6 presents the evolution of the deformation of the solid-fluid interface in the time interval t ∈ [0, 500] s, in a similar manner as in References 21,50. The red curve in Figure 6 indicates the equilibrium deformation of the solid-fluid interface obtained in test case 1; cf. Figure 4. The surface tension carried by the fluid-fluid interface exerts a localized load on the solid-fluid interface at the contact line. In the sharp-interface limit → +0, this concentrated load is locally only balanced by surface tension, not by elasticity 21 or inertia. Hence, the contact-line traction leads to instantaneous formation of a kink 51 in the surface of the soft substrate. Similar behavior is observed in Figure 6 for the diffuse-interface model, by virtue of the fact that is significantly smaller than the elasto-capillary length. Due to the finite thickness of the diffuse fluid-fluid interface, the kink is, however, obtuse and exhibits finite curvature. The kink develops into a wetting ridge of approximately 8.75 m height; see also Figure 4.
To further elucidate the dynamics of the solid-fluid interface, Figure 7 displays the evolution of the maximum displacement, the minimum displacement and the displacement at the center of the substrate. The maximum displacement corresponds to the position of the apex of the wetting ridge. The apex height increases monotonously. After 800 s, the wetting ridge has reached a height of 8 m, corresponding to more than 90% of its equilibrium height. The minimal displacement decays monotonously. After 800 s, the deviation of the minimum displacement relative to its equilibrium value is still significant. The center displacement displays non-monotonous behavior. Initially, the indentation below the

4.3
Convergence of the adaptive-refinement procedure Finally, we assess the a posteriori error-estimation procedure, and the correspond adaptive-refinement process outlined in Section 3.3. For this purpose, we regard the sessile-droplet test case with parameter set #3 in Table 1. The parameters are identical to parameter set #2, except for the refinement fractionsû,̂and̂, which have been set to 0.9,0.9 and 0.5 instead of 0,0.5 and 0, respectively. The reason for modifying these refinement fractions is that targeted refinements based on the fluid-velocity field and the chemical potential do not have a visible effect on the results, but do have an effect on the (asymptotic) convergence rates of the adaptive method. The precise values of the refinement fractions have been selected heuristically. Figure 8 plots the error estimates for the fields f ∈ {û,p,̂,̂,̂,̂} versus the dimension of the approximation space for the respective fields, #dof f L , for refinement levels L = 0, 1, 2, … , L max − 1. Note that the two-level hierarchical F I G U R E 8 Convergence behavior of the adaptive-refinement procedure: Error estimates for the fields f ∈ {û,p,̂,̂,̂,̂} vs the dimension of the adaptively refined approximation space (solid) and upper bound according to (39) (dashed). Triangles indicate first approximation after -continuation process error estimate is not available for the ultimate refinement step, L max . We present results for time step n = 2, but the convergence properties of the adaptive method are essentially independent of the time step. The dashed lines indicate the upper bounds for the respective fields according to the inner summand in (39). The triangles mark the Kth refinement step, that is, the first approximation with the prescribed interface thickness after the initial -continuation process; see Section 3.4. Figure 8 conveys that the error estimates for all fields decay monotonously. The error bounds, composed of the support-wise error indicators according to (38), are essentially proportional to the error estimates. This corroborates the suitability of the error indicators for guiding the adaptive refinements. The effectivity indices of the error bounds, that is, the ratios of the error bounds to the actual error estimates, are consistently below 10. This is appropriate, given that many error contributions are duplicated in the error bound on account of the overlap of the basis-function supports. Figure 8 does not provide an unambiguous asymptotic convergence rate for the various fields, especially for the deformation and the solid pressure. The approximations of the binary-fluid fields,û,p,̂, and̂, appear to convergence as #dof −3 . In general, the optimal convergence rate for the approximation of a d-variate function u by a piece-wise polynomial approximation, u h,k , of degree k is: The approximations ofû,̂, and̂are piece-wise cubic (k = 3), while the Sobolev index for these fields is m = 1. The approximation of the fluid pressure,p, is piece-wise quadratic (k = 2), and the corresponding Sobolev index is m = 0. For all the binary-fluid fields, the optimal convergence rate according to (47) is therefore −3∕2, instead of the observed rate −3. According to (47), the observed rate of convergence as provided by the adaptive algorithm for the binary-fluid fields apparently pertains to a one-dimensional object. This is consistent with the fact that, essentially, the adaptive algorithm can very effectively reduce the approximation errors by introducing refinements only on the diffuse interface. Because the diffuse-interface appears as a manifold of co-dimension 1 at coarse levels of refinement, the convergence behavior of the adaptive procedure resembles the optimal convergence rate in 1 dimension. We conjecture that at even higher levels of refinement, the asymptotic convergence rate according to (47) is recovered. The results in Figure 8 are inconclusive with respect to the convergence rate of the adaptive algorithm for the solid fields, including the deformation.

CONCLUSION
We presented an adaptive IGA methodology for simulating the elasto-capillary FSI of a binary fluid with an elastic solid, based on an incompressible NSCH diffuse-interface model for the binary fluid, a hyperelastic model for the solid, and an ALE formulation for the FSI problem. We considered numerical approximations constructed from adaptive high-regularity THB splines, as employed in the IGA paradigm. The adaptive strategy relies on a two-level hierarchical a posteriori error estimate. The hierarchical error estimate directs a support-based refinement procedure, which naturally accommodates the larger basis-function supports of the high-regularity THB splines relative to traditional C 0 -continuous finite-element approximations. To avoid the nonrobustness of partitioned solution methods for elasto-capillary FSI, we applied a fully monolithic solution procedure to solve the aggregated BFSI problems in each time step and each iteration of the adaptive-refinement process. The required derivatives were obtained by means of automatic differentiation. To further enhance the robustness of the monolithic solution procedure, we introduced a so-called -continuation approach, in which the value of the interface-thickness parameter, , is artificially enlarged and reduced to its prescribed value in the initial iterations of the adaptive procedure.
To illustrate the properties of the adaptive IGA methodology for elasto-capillary FSI, we presented numerical computations for a configuration pertaining to a sessile droplet on a soft solid substrate. Comparison of the computed results with experimental data for the considered test case exhibited excellent agreement. In terms of agreement with the experimental data, the present results provide a significant improvement relative to previous results for the same test case, by virtue of modified soft-solid constitutive behavior, namely neo-Hookean in the present work versus St. Venant-Kirchhoff in previous studies. The overall good agreement between the computed results and experimental data corroborates the potential of computational BFSI models based on the NSCH equations for describing elasto-capillary FSI phenomena.
We observed that the error-estimate-based adaptive-refinement procedure automatically concentrates the refinements in the vicinity of the diffuse interface and near the apex of the wetting ridge. By virtue of the excellent approximation properties of the considered high-regularity THB splines, the diffuse interface is essentially resolved with only a few knot spans across the diffuse interface. For the considered test cases, the adaptive-refinement method, guided by the two-level hierarchical error estimate, provided very good monotonous convergence behavior. We observed super-convergent behavior of the adaptive method for the approximations of the binary-fluid fields, which is conjecturally due to the fact that the diffuse interface appears as a manifold of co-dimension one.
We believe that the presented computational framework can be further developed by considering alternate error-estimation methods. The two-level hierarchical error estimate provides a reliable and robust estimate of the error, but carries high computational complexity. Moreover, on account of the many disparate fields involved in BFSI, the selection of refinement fractions, that is, the relative amount of refinement to be assigned to each field, is nonobvious. In addition, the robustness of the computational framework requires further enhancement for high viscosity-ratio binary fluids or, in particular, binary fluids with low-viscosity components, for example, water-air mixtures. Several measures have been introduced in this work to improve the robustness of the methodology, namely, soft-clipping of the mixture-density interpolation, the Arrhenius constitutive relation for the mixture viscosity, a skew-symmetric reformulation of the convective trilinear form, a fully monolithic solution procedure and -continuation, but further enhancements are necessary.