A 3D isogeometric BE–FE analysis with dynamic remeshing for the simulation of a deformable particle in shear flows

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. 326 (2017) 70–101 www.elsevier.com/locate/cma

A 3D isogeometric BE–FE analysis with dynamic remeshing for the simulation of a deformable particle in shear flows Jorge Maestre a , Jordi Pallares a , ∗, Ildefonso Cuesta a , Michael A. Scott b a Departament d’Enginyeria Mecànica, Universitat Rovira i Virgili, Av. Països Catalans, 26, 43007 Tarragona, Spain b Department of Civil and Environmental Engineering, Brigham Young University, Provo, UT 84602, USA

Received 31 March 2017; received in revised form 3 August 2017; accepted 4 August 2017 Available online 16 August 2017

Highlights • • • • •

A T-spline-based isogeometric approach for demormable capsules is presented. An IGA-BIE is derived for the computation of the flow in pipes of arbitrary sections. An adaptive refinement scheme is formulated to handle the mesh quality on the capsule. The accuracy and stability of the method are checked with some benchmark examples. The behavior of a capsule flowing through a real geometry of a capillary is simulated.

Abstract A three-dimensional isogeometric coupled boundary element and finite element approach based on analysis suitable T-splines is developed for the simulation of deformable capsules suspended in shear flows. Boundary element analysis is used to solve the fluid Stokes equation whereas the hydrodynamic membrane load is computed via isogeometric analysis under the assumption that the membrane is a hyper-elastic thin shell with negligible bending resistance. The smoothness of the T-spline basis functions accommodate large deformations of the capsule without the need for additional smoothing techniques, and can be used to accurately compute the membrane load. A balanced distribution of membrane elements can be constructed using an unstructured locally refined mesh. These properties are coupled with an adaptive temporal implicit integration scheme. Several benchmark examples are solved to illustrate the accuracy and potential of the method. The approach is then applied to simulate the dynamics of a capsule in a real geometry of a brain capillary. c 2017 Elsevier B.V. All rights reserved. ⃝

Keywords: Isogeometric analysis; BEM; FEM; T-spline; Deformable capsule; Dynamic refinement

∗ Corresponding author.

E-mail address: [email protected] (J. Pallares). http://dx.doi.org/10.1016/j.cma.2017.08.003 c 2017 Elsevier B.V. All rights reserved. 0045-7825/⃝

J. Maestre et al. / Comput. Methods Appl. Mech. Engrg. 326 (2017) 70–101

71

1. Introduction The hydrodynamic movement of suspended capsules in a flowing fluid is an important area of research which can be applied to a variety of industrial and biological applications. In the biomedical area, the study of the dynamic behavior of deformable vesicles and Red Blood Cells (RBC) has led to a better understanding of thrombosis, thalassemia and other microcirculatory diseases [1]. In pharmacological treatments, artificial capsules are used as an effective technique for drug delivery. Additionally, micro-capsules are being applied in cosmetics, agriculture and the food industry [2]. The characterization and study of deformable capsules is required for the efficient design and optimization of different processes connected with the capsule dynamics. Experiments are indispensable but, because of the small length scales involved, they are difficult and costly. Numerical models and methods have opened new perspectives and they offer a powerful tool to analyze capsule deformation. A capsule consists of a viscous internal fluid enclosed and protected by a flexible membrane. The interaction between the capsule and the external flowing fluid is complex. Analytical solutions have been proposed for deformable capsules in simple shear flows [3,4], but the application is limited to small deformations. Capsules usually experience large deformations which change the shape and modify the behavior of the flowing fluid. In these cases, the only option is a numerical model. Different models have been successfully used to solve the fluid–structure interaction. The majority of these models are based on Finite Element Methods (FEM) [5,6], Finite Difference Methods (FDM) [7], Immersed Boundary Methods (IBM) [8–12] and Lattice Boltzmann Methods (LBM) [13–16]. These approaches are relatively well established and are applicable to a wide range of problems (i.e. different geometries), Reynolds numbers and fluid properties. In many microfluidic applications, the Reynolds number is very low and the inertial effects are negligible. Assuming a Newtonian or nearly-Newtonian fluid, the flow is governed by the Stokes equation. In this case, the Boundary Element Method (BEM) is a favorable alternative for several reasons: (1) only surface meshes are needed; (2) Computer Aid Design (CAD) surfaces can be accommodated directly using isogeometric analysis (IGA); and (3) the Sommerfeld condition is implicitly satisfied in the formulation, which simplifies setting the boundary conditions in unbounded flows (the reader is referred to [17] for more details about the Sommerfeld condition). Several BEM models for deformable capsules have been reported in the literature. Pozrikidis [18] studied the transient large deformation of a three-dimensional spherical isoviscous capsule in a shear flow using a global geometric representation of the membrane based on quadrilateral elements. Ramanujan and Pozrikidis [19] extended the formulation to unstructured meshes built with quadratic triangular elements and then investigated the large deformation of various capsule shapes considering the effect of different fluid viscosities. Leyrat-Maurin and Barth`eesBiesel [20] and Qu´eguiner and Barth`es-Biesel [21] developed an axisymmetric BEM for the analysis of deformable capsules through pores, tubes and constrictions. Diaz et al. [22], Diaz and Barth`e-Biesel [23] and Lefebvre and Barth`es-Biesel [24] improved the previous formulation by using cubic B-splines and then they analyzed the influence of various parameters in the deformation of the capsules (such as membrane properties and shapes, fluid viscosity, and membrane pre-stress). Lac et al. [25] and Lac and Barth`es-Biesel [26] presented a three dimensional formulation based on bicubic B-spline basis functions for deformable capsules in a planar shear and a hyperbolic flow. This method achieves a high-order of approximation and continuity of the physical variables leading to accurate results. Deformable capsules in moderate and high shear flow were also tackled by Dodson and Dimitrakopoulos [27,28] who developed a spectral BEM. More recently, accelerated BEM has been developed by Zhu et al. [29] to study the motion of a capsule in a wall-bounded oscillating shear flow and by Touchard et al. [30] to study the deformation of a capsule flowing through a micro-channel with a localized constraint. Following the pioneering works of Zarda et al. [31] and Skalak et al. [32,33], most of studies consider the membrane to be a thin hyper-elastic shell with negligible bending resistance. This approach is applicable to a wide range of artificial capsules and biological cells and it is in agreement with experimental findings that confirm the relatively low bending resistance of the membrane [34–36]. It is worth mentioning that other more detailed models have been developed which take into account the bending resistance, surface incompressibility, membrane viscosity and rheology (see for example [37–41]), although they are computationally intensive. A complete review of the mechanical properties of deformable capsules can be found in [42]. The accuracy of the BEM depends on the computation of the hydrodynamic load over the membrane. Using local equilibrium, a high order differentiation of the geometry is required. In a traditional BEM, the basis functions are C 0 continuous between element, thus the spatial derivatives are discontinuous across them. Pozrikidis [18] and Ramanujan and Pozrikidis [19] found that instabilities associated with C 0 basis functions negatively affect to the

72

J. Maestre et al. / Comput. Methods Appl. Mech. Engrg. 326 (2017) 70–101

solution. They then applied a smoothing scheme to relocate the nodes and an ad-hoc averaging technique of the membrane stress to alleviate these numerical instabilities. A smoothing technique was also employed by Dodson and Dimitrakopoulos [28] to force the continuity of the first derivatives between spectral elements. Lac et al. [25] reported instabilities near geometric poles despite the high order continuity of the B-spline functions. To overcome this issue, Walter et al. [43] formulated a coupled BEM–FEM model for deformable capsules. The loads are computed through a variational principle thus reducing the continuity requirements of the basis functions and producing more stable results. This formulation has been recently used to study the dynamics of deformable capsules [44–48]. IGA was introduced by Hughes et al. [49] as a generalization of FEM. In IGA, the basis functions, which define the geometry, are used as the basis for analysis thus establishing a tight link between the CAD geometry and the engineering analysis. IGA possesses several distinct advantages: (1) Exact geometry is embedded in the numerical analysis; (2) High order approximation and continuity of the bases; (3) Generalized refinement schemes which exactly preserve the geometry. Most initial studies were based on Non-uniform rational B-splines (NURBS). NURBS-based IGA has experienced widespread acceptance in the engineering analysis community. In biological applications of IGA, several studies have demonstrated superior accuracy and robustness when compared to finite element analysis (FEA) (see [50–53]). The extension of IGA to BEM was initiated by Politis et al. [54] and Simpson et al. [55], who applied the technique to potential and elasto-static problems. In the last few years, this method has grown rapidly [56–61], although the application to biomedical flows and deformable capsules is still in its infancy. Joneidi et al. [62] presented an IGA-BEM model to study the dynamic behavior of drops and inextensible capsules in a shear flow. Heltai et al. [63] formulated a non singular IGA-BEM for steady three dimensional Stokes flows with the idea of extending model to flows induced by biological organisms as swimming bacteria. Recently, Heltai et al. [64] presented a coupled IGA BEM and shell approach for the analysis of the interaction between thin elastic structures and Stokes flows. Unfortunately, NURBS-based IGA has several important limitations. NURBS are restricted to rectangular parametric domains. This greatly limits their use when building complex geometric models. In the case of spherical surfaces, the presence of degenerate elements around poles is common and can introduce numerical problems in IGA. T-splines appeared as a more advanced technology to overcome these limitations [65]. Unlike NURBS, unstructured meshes can be naturally accommodated and the resulting surface is locally refinable. Scott et al. [66] and Li et al. [67] developed analysis suitable T-spline which ensures the appropriate mathematical properties for IGA. Borden et al. [68] and Scott et al. [69] introduced the concept of B´ezier extraction as an efficient and elegant technique to incorporate NURBS and T-spline bases into FEA and Scott et al. [70] extended the analysis suitable T-spline to BEM. The T-spline-based IGA analysis has been successful applied to a variety of fields such as contact [71], acoustics [72], gravity waves [73,74], elasticity [75,76], fracture and damage [77–81], fluid– structure interaction [82,83], electromagnetic [84] and tumor growth [85]. We should mention that alternatives to NURBS and T-splines have been introduced such as PHT-splines [86], hierarchical splines (B-splines, NURBS and T-splines) [76,87,88], truncation splines (B-spline and NURBS) [89,90] and Truncated Hierarchical Catmull–Clark Subdivision [91], among others. The objective of this study is to apply IGA, based on the analysis suitable T-spline, to study the dynamics of deformable capsules in suspended flows. The BEM is used to solve the fluid dynamics problem which is coupled to a structural membrane. The membrane is considered to be a thin elastic shell, thus the discretization of the interface is shared for both methods. The smoothness of the basis prevents the saw-tooth effect (i.e. the artificial sharp edge formation across elements, see [74]) and makes the method stable for large deformations of the capsule without the need for artificial smoothing techniques. Degenerate elements at the poles can be avoided by using unstructured T-splines. In addition, invoking the local refinability of the T-spline, the mesh can be dynamically adapted to the specific problem. These properties give the analysis suitable T-spline-based IGA the potential to be an efficient, accurate and robust tool for the analysis of deformable capsules. This work is structured as follows. In Section 2 we show the governing equations of the fluids and of the elastic membrane. The numerical aspects of IGA applied to deformable capsules is discussed in Section 3. In Section 4 various numerical examples are shown to reveal the potential of the approach. Finally, we outline the conclusions of this work in Section 5. 2. Formulation of the problem 2.1. Hydrodynamics We consider a capsule defined by the boundary ΓC and confined in a tube ΓT = Γ I ∪ Γ O ∪ ΓW , where Γ I , Γ O and ΓW are the inlet, outlet and wall boundaries, respectively. The tube domain, Ω1 , is filled with an external fluid.

J. Maestre et al. / Comput. Methods Appl. Mech. Engrg. 326 (2017) 70–101

73

Fig. 1. Sketch of the problem.

The capsule is formed by a thin deformable membrane whose internal domain, Ω2 , contains an internal fluid. The capsule moves freely by the action of the external flow in the tube, while the internal fluid moves in response to the deformation of the membrane. Both fluids are assumed to be incompressible and Newtonian, with the same density ρ and different viscosities, µ1 = µ and µ2 = λµ, where λ is a viscosity ratio. A sketch of the problem is shown in Fig. 1. In microcirculation problems the Reynolds number is small, Re ≪ 1, and the viscous effects are assumed to be dominant compared to the inertial effects. In these cases, the internal and external fluids can be described by the Stokes equations, { ∇ · vi = 0 in Ωi , i = 1, 2 (1) ∇ · σ i = −∇ pi + µi ∇ 2 vi = 0 where σ i is the fluid stress tensor given in terms of the velocity vi and pressure pi in the fluid domain Ωi . In this work we consider only non slip conditions. This ensures continuity of the velocity across the membrane. The stress on both sides of the membrane is allowed to be discontinuous. The difference in the stress across the membrane is the stress jump, ∆f = (σ 1 − σ 2 )n, where n is the unit normal vector pointing outwards of the capsule. The stress jump depends on the mechanical properties of the membrane as described in the next section. 2.2. Mechanical properties of the membrane The dynamical and mechanical behavior of capsules is governed by the micro-mechanical properties of the membrane and the properties of the external and internal fluids. Biological membranes are complex systems that can experience large deformations, which change the shape and the mechanical characteristics of the capsule. Various models at different scales and levels of detail have been proposed depending on the type of biological capsule and the application [40,92–99]. In this study, we employ a general continuous elastic model. As proposed by several authors (see for example [22,24,25,28,31,43,94]), we assume that the membrane is extremely thin and the strain and stress are constant through the thickness of the membrane. Therefore, the membrane behaves like a bi-dimensional elastic sheet and the bending resistance is neglected. We describe the membrane by a curvilinear coordinate system defined by two parametric coordinates ξ = (ξ 1 , ξ 2 ). Thus, the reference configuration of the membrane, i.e. t = 0, is given by X(ξ ). The membrane in the current configuration is denoted by x(ξ , t). In other words, uppercase and lowercase letters correspond to the reference and current configuration, respectively. In the current configuration, the local covariant bases are constructed as follows: g1 × g2 g1 = x,1 , g1 = x,2 , g3 = n = (2) |g1 × g2 | ∂(·) where the subscript (·),α = ∂ξ α denotes a partial derivative with respect to the parametric coordinates (see Fig. 2). The 1 2 3 local contravariant bases (g , g , g ) are defined such that gi · g j = δ ij , where δ ij is the mixed Kronecker delta. The metric tensors associated with each basis are gi j = gi · g j and g i j = gi · g j . The local covariant, contravariant and metric tensor in the reference configuration are analogously defined and represented by uppercase letters. The local equilibrium equation of the membrane in the current configuration is given by

∇s · τ + t = 0

(3)

74

J. Maestre et al. / Comput. Methods Appl. Mech. Engrg. 326 (2017) 70–101

3 ). Fig. 2. Curvilinear coordinate system referred to the membrane (covariant bases {gi }i=1

where we assume that the inertia of the membrane is negligible, which is reasonable in the case of thin membranes. τ is the in-plane Cauchy stress tensor that corresponds to the force per unit arc length of the membrane, ∇s = (I−n⊗n)·∇ is the surface gradient operator, being I the identity matrix, and t = ∆f is the traction on the membrane that is equal to the jump stress in Eq. (3). We assume that the membrane is a hyper-elastic isotropic material. This material is a subclass of a general elastic material in which the constitutive equation can be derived from a strain energy function Ψ (F): τ=

∂Ψ (F) FT ∂F Js

(4)

where the tensor F = dx/dX = gi ⊗ Gi is the deformation gradient (Gi and Gi are the contravariant and covariant bases in the reference configuration, respectively) and Js = |F| is the surface dilatation. In the case of isotropic bi-dimensional models, the deformation gradient and the strain energy function only depend on two strain invariants: I1 = G αβ gαβ and I2 = |G αβ ||gαβ | − 1 = Js2 − 1

(5)

Note that the Greek indices α and β go from 1 to 2, while Latin indices i and j go from 1 to 3. The Cauchy stress tensor in the contravariant form can then be written as ∂Ψ (I1 , I2 ) αβ ∂Ψ (I1 , I2 ) αβ τ αβ = 2Js−1 G + 2Js g (6) ∂ I1 ∂ I2 In the literature, several constitutive laws have been proposed to take into account the mechanical properties of biological membranes [100]. In this work, two common models are adopted: the Neo-Hookean model and the model proposed by Skalak et al. [32]. The bi-dimensional Neo-Hookean model is usually used to describe the behavior of proteinreticulated membranes. The energy function is given by ) ( G sN H 1 NH −1 (7) Ψ = I1 + 2 I2 + 1 where G sN H is the shear modulus. The model reported by Skalak et al. [32] is appropriate to represent the behavior of red blood cell membranes. The energy function is given by ) G sSk ( 2 I1 + 2I1 − 2I2 + C I22 (8) 4 in which G sSk is the shear modulus and C > −1/2 is a parameter that measures the ratio between the area dilatation and the shear modulus such that K sSk = G sSk (1 + 2C). In contrast to the softening Neo-Hookean model, the Skalak law shows strain-hardening behavior and for C ≫ 1 the model is almost area-incompressible, as observed in RBCs. Ψ Sk =

3. Formulation of the method 3.1. Surface discretization The numerical model is defined by the external boundary surface ΓT and the internal boundary surface ΓC that represent the tube and the capsule, respectively. In this work, both surfaces are discretized using cubic T-splines in the framework of isogeometric analysis. T-splines provide a high order surface parametrization and allow local refinements. The continuity and smoothness of the surface can be directly controlled via knot vectors, even during

J. Maestre et al. / Comput. Methods Appl. Mech. Engrg. 326 (2017) 70–101

(a) T-mesh.

75

(b) Physical space.

Fig. 3. (a) T-mesh construction of a (b) T-spline surface. The red lines delimit the T-mesh elements in the T-mesh, being the dashed lines, lines of reduced continuity. In (a) the dark highlighted quadrilateral denotes the T-mesh element number 1 and it is indicated in (b) in the physical space (soft highlighted quadrilateral). The highlighted circles in (a) correspond to non-zero basis functions over the T-spline element 1. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

the refinement process. In addition and contrary to the traditional C 0 FEA, the use of T-splines allows the unique evaluation of the normal and the curvature on smooth surfaces. These features lead to higher order accuracy and stability in dynamic problems. T-splines are a generalization of NURBS that allow for hanging-nodes and extraordinary points. The surface is constructed form an unstructured mesh with local refinements. We briefly introduce T-spline concepts. The interested reader is referred to [65,66,75,101] for additional details. A T-spline surface is defined by a polygonal mesh called T-mesh (see Fig. 3(b)). The T-mesh is constructed from a tiling of quadrilateral elements. Using the T-mesh, the T-spline basis functions can be inferred by constructing a set of local knot vectors over the mesh. A control point P A ∈ ℜ3 and weight w A ∈ ℜ, A = 1, 2, . . . , m cp , is associated with each basis function, where m cp is the number of control points. The T-spline surface that maps the parametric ˆ → Γ , can be expressed as space to the physical space, x : Γ ∑m cp ˆ A=1 P A w A N A (ξ A ) ˆA x= ∑ ξˆ A ∈ Γ (9) m cp ˆ A) w N ( ξ A A A=1 ˆA . where N A (ξˆ A ) ∈ ℜ≥0 is a T-spline basis function restricted to the parametric domain Γ In isogeometric analysis it is common to describe a T-spline from a finite element point of view [69]. The physical ˜. Over each element a small number domain is divided in a series of elements, xe (ξ˜ ), defined over a parent domain Γ of basis functions are non-zero (see Fig. 3(b)). We construct a correspondence array, A = I E N (a, e), that maps a local basis function index a in element e to a global basis function index A. The T-spline representation can then be written as (Pe )T We Ne (ξ˜ ) ˜ xe (ξ˜ ) = ξ˜ ∈ Γ (10) (we )T Ne (ξ˜ ) me me where Pe = {Pae }a=1 is the vector of control points, we = {Wae }a=1 is the vector of weights, We = diag(we ) is the m e diagonal matrix, Ne (ξ˜ ) = {Nae (ξ˜ )}a=1 is the vector of basis functions supported by the element and m e is the number of non-zero basis functions defined over the element e. Note that the number of non-zero basis functions can change from element to element. Furthermore, a basis function is supported by several elements and its properties depend on

76

J. Maestre et al. / Comput. Methods Appl. Mech. Engrg. 326 (2017) 70–101

Fig. 4. Representation of (a) a structured mesh with degenerate elements and (b) a unstructured mesh with an extraordinary point. The pole is marked with a red circle and the extraordinary point with a red triangle. Each color strip corresponds to one ring of neighbor elements.

the topology of the T-spline mesh. As proposed by Scott et al. [69], we use the B´ezier extraction concept as an elegant method which standardizes the T-spline framework for isogeometric analysis. B´ezier extraction generates a linear operator that maps the T-spline basis over each element to a canonical set of mb Bernstein polynomials. For the cubic case, we define a set of bivariate cubic Bernstein polynomials, B = {Bb (ξ˜ )}b=1 , where m b = 16. The T-spline basis can then be expressed as Ne = Ce B

(11)

where Ce ∈ ℜm e ×m b is called the matrix extraction operator corresponding to element e. Introducing Eq. (11) in Eq. (10), each element can be written in terms of the canonical set of Bernstein polynomials as, xe (ξ˜ ) =

(Qe )T Wb,e B(ξ˜ ) (Pe )T We Ce B(ξ˜ ) = (we )T Ce B(ξ˜ ) (wb,e )T B(ξ˜ ) m

˜ ξ˜ ∈ Γ

(12) m

b b where Qe = {Qeb }b=1 is a vector of B´ezier control points and wb,e = {Wbe }b=1 is a vector of B´ezier weights. Eq. (12) can be expressed in reduced form as

xe (ξ˜ ) = (Pe )T Re (ξ˜ ) = (Re )T (ξ˜ )Pe

˜ ξ˜ ∈ Γ

(13)

is a vector of rational basis functions in B´ezier form. where R (ξ˜ ) = In this paper we deal with surfaces that have extraordinary points and degenerate elements (see Fig. 4). In both cases, the continuity of the geometry is reduced locally. An extraordinary point is a vertex which is adjacent to a number of edges different from four. In the present work, we leverage the constrained B´ezier extraction framework reported in [70] around extraordinary points. This approach enforces G 1 continuity in the one-ring neighborhood of elements adjacent to an extraordinary point. A degenerate element is an element in which one of the edges of the quadrilateral element is collapsed to a point forming a triangular element. There may be several elements that collapse to the same point. We will call these degenerate points poles (e.g. the poles of a sphere). The continuity across a pole can be chosen to be C 0 if a continuous basis function is constructed as the sum of all the coincident basis functions. However, it is not possible to compute the surface normal and curvature at the poles and this can affect the accuracy and robustness of numerical results [72,102,103]. To overcome this issue, we perform a simple transformation to obtain non-singular elements. The procedure is as follows (see Appendix for the complete matrix expressions): e

me {Rae }a=1

1. Given a degenerate cubic quadrilateral B´ezier element defined by the homogeneous coordinates Qw = mb {wib Qi , wib }i=1 , the element is approximated by a triangular cubic B´ezier element Tw . The transformation can be expressed in a matrix form as Tw = C QT Qw

(14)

where C QT ∈ ℜm T ×m b transform quadrilateral B´ezier control points to triangular B´ezier control points and m T = 9 is the number of control points of the cubic triangular B´ezier element. Note that this step only

J. Maestre et al. / Comput. Methods Appl. Mech. Engrg. 326 (2017) 70–101

77

approximates the original T-spline surface. We can apply a subdivision scheme to manage the tolerance of this approximation [104]. ˆw , Q ˆw , Q ˆw ) 2. The cubic triangular B´ezier element is then divided into three cubic quadrilateral B´ezier elements (Q 1 2 3 such that ˆiw = CT Q Tw Q i

i = 1, 2, 3

(15)

where CiT Q ∈ ℜm b ×m T converts from triangular B´ezier patches to quadrilateral B´ezier patches. This procedure was originally developed by the CAGD community [104,105]. In the context of B´ezier-extraction, this procedure is easy to implement since it consists of simple matrix multiplications of extraction operators with the conversion matrices. 3.2. Isogeometric boundary element The boundary integral formulation for Stokes flows is well known [106]. We focus on those aspects related to isogeometric discretizations. In particular we use a collocation approach as described in [70]. The field velocity and fluid stress are discretizated in a finite dimensional space V h defined by the T-spline basis functions. Using the global indexing of basis functions, both fields can be expressed as v(x) =

m cp ∑

R A (x)v A

A=1 m cp

f(x) =

∑

(16) R A (x)f A

A=1 m

m

cp cp where {v A } A=1 and {f A } A=1 are the velocity and fluid stress control vectors corresponding to the global T-spline basis m cp {R A } A=1 . Given a continuous surface Γ , the single and double-layer potential operators at the collocation point x0 are defined as ∫ ∫ 1 1 − Ti j (x, x0 )v j (x)dΓ (x), (G Γ f)(x0 ) = G i j (x, x0 ) f j (x)dΓ (x) and (HΓ v)(x0 ) = (17) 8π µ Γ 8π Γ

respectively. In discrete form, these operators can be written as m cp ∫ 1 ∑ (G Γ f)h (x0 ) = G(x, x0 )R(x) A f A dΓ 8π µ A=1 Γ and (18) m cp ∫ ∑ 1 − T(x, x0 )R A (x)v A dΓ (HΓ v)h (x0 ) = 8π A=1 Γ ∫ Although both integrals are well defined, the symbol − indicates that the double-layer integral is evaluated in the Cauchy Principal Value sense (CPV) when the collocation point lies on the integration domain Γ . The tensors G(x, x0 ) and T(x, x0 ) are the free-space Green’s functions for Stokes flow. In three dimensions, these tensors are given by: G i j (x, x0 ) Ti j (x, x0 )

ri r j δi j + 3 |r| |r| ri r j r k n k = −6 |r|5 =

i, j and k = 1, 2, 3

(19)

where δi j is the Kronecker delta and r = x − x0 . We view the Stokes problem as a superposition of two flows. In other words, the total fluid velocity and fluid traction vectors can be expressed as v = v∞ + v D and f = f∞ + f D , respectively. The unperturbed flow, denoted by the superscript (·)∞ , corresponds to a known pipe flow without the capsule (see Fig. 1). Whereas the disturbed flow,

78

J. Maestre et al. / Comput. Methods Appl. Mech. Engrg. 326 (2017) 70–101

denoted by the superscript (·) D , is caused only by the presence of the capsule. The boundary integral equation (BIE) applied to the capsule can be expressed as [I − (1 − λ)C(x0 )] v(x0 ) = (1 − λ)(HΓC v)(x0 ) − (G ΓC ∆f)(x0 ) D

∞

− (G ΓT f )(x0 ) + v (x0 )

(20)

x0 ∈ ΓC

and in the tube as 0 = (1 − λ)(HΓC v)(x0 ) − (G ΓC ∆f)(x0 ) − (G ΓT f D )(x0 ) x0 ∈ ΓT

(21)

where I ∈ ℜ3×3 is the identity matrix and the term C(x0 ) ∈ ℜ3×3 is the solid angle which can be computed directly using the formulation proposed by Mantic [107] or indirectly appealing to the “rigid body motion” method [108]. We employ the rigid body motion technique since it leads to a regularization of the second-layer integral and avoids the computation of the integral in the CPV sense. Thus, the second-layer potential term on the capsule can be written as ∫ 1 ˆ T(x, x0 )(v(x) − v(x0 ))dΓ (x) (22) (HΓC v)(x0 ) = C(x0 )v(x0 ) + (HΓC v)(x0 ) = 8π ΓC which in discrete form can be written as m cp ∫ 1 ∑ h ˆ (HΓC v) (x0 ) = T(x, x0 ) (R A (x) − R A (x0 )) v A dΓ 8π A=1 Γ

(23)

The single-layer and regularized double-layer potential operators are then weakly singular. As proposed by Lachat and Watson [109], we address this singularity by making a local change of variable in the element that contains the collocation point. In Eqs. (20) and (21) we assume that the inlet Γ I and outlet Γ O boundaries (see Fig. 1) are sufficiently far from the capsule such that the perturbed velocity flow is practically negligible. The disturbed fluid stress on the inlet and outlet can be retained as unknown variables in the BIE or explicitly computed by an auxiliary equation as proposed by Pozrikidis [95]. This approach reduces the number of degrees of freedom (DoF) of the BIE and it is used in this paper. Approximating the disturbed fluid stress at the inlet as fΓDI ≃ −∆ p D n and setting fΓDO = 0 at the outlet, the disturbed pressure drop, ∆ p D , can be calculated as ∫ 1 D ∆p = ∞ ∆f(x) · v∞ (x) + (1 − λ)f∞ (x) · v(x)dΓ (x) (24) Q ΓC where Q ∞ ≃ Q is the total flow rate. We use generalized Greville abscissae as the collocation points of the BIE [70]. This strategy accommodates extraordinary points and hanging-nodes (i.e. T-junctions) and it has been shown that it provides suitable results in the context of isogeometric analysis. The collocated BIE can be written as [ ][ ] [ ][ ] ˆCC − RT (xC ) 0 vC GCC GC W ∆f (1 − λ)H = D 0 GW C GW W fW (1 − λ)HW C 0 (25) [ ] [ ∞ ] GnC v (x ) − ∆ pD − ∞ C GnW v (xW ) where H and G denote the single and double-layer global matrices, respectively. The first subscript refers to the collocation points (on the capsule, (·)C , or on the tube wall, (·)W ) and the second subscript refers to the integration domain. Gn is a global vector that denotes the single-layer operator applied to the inlet cap assuming the approximation of a constant fluid stress vector in the normal direction. R is a matrix that contains the rational basis functions evaluated at the collocation points, xC . The vector v∞ is the unperturbed velocity vector evaluated at the collocation points. Note that the jump stress vector on the capsule, ∆f, can be determined through the membrane equilibrium equation (Eq. (3)). The computation of the membrane stress will be explained in the next section. The D velocity vector of the capsule, vC , the disturbed fluid stress vector of the tube wall, fW , and disturbed pressure drop, D ∆ p , are unknown variables. The problem is completely defined by adding the pressure drop equation (Eq. (24)). Using standard approaches, the global matrix and vectors can be constructed by integrating and assembling each B´ezier element into the global system [69]. We can distinguish three different kinds of integrals: not singular (or regular), nearly-singular and singular, if the collocation point is placed outside, close or inside of the integration

79

J. Maestre et al. / Comput. Methods Appl. Mech. Engrg. 326 (2017) 70–101

element, respectively. Since singular integrals are properly regularized and treated, the integrals can be numerically evaluated by a quadrature rule. As commonly used in the traditional BEM [108,110], we employ a Gauss–Legendre quadrature rule. This quadrature has been also extended satisfactorily in IGA (see for example [63,70,72,102]). However, it must be noted that the Gauss–Legendre quadrature exhibits a sub-optimal converge, as shown for example in [63,111], and it does not take into account the high continuity of the T-spline bases across the elements. New schemes have been proposed in the context of IGA (see [111,112]), but they require more complex algorithms and they are mesh-depending. The evaluation of the nearly-singular integral is also a challenge and it is an active research area. Several methods have been presented on that subject based on adaptive integration schemes [102], regular element subdivision [109], analytical evaluation of the singular part [113] and transformation techniques [114,115]. In this paper, we use a relatively large number of Gauss points to address these issues and to obtain good accuracy. The reader is referred to Maestre et al. [74] for more details about the numerical evaluation of these integrals. 3.3. Membrane discretization As proposed by Walter et al. [43], we leverage a weak formulation of the membrane mechanics. This technique, as opposed to the collocation technique, relaxes the continuity requirements of the basis functions and gives accurate and stable results without the need to implement ad-hoc smoothing techniques. We present this formulation in the context of isogeometric analysis. The principle of virtual work (PVW) applied on the membrane in the current configuration can expressed as ∫ ∫ δu(x) · t(x)dΓ (x) = δe(x) : τ (x)dΓ (x) (26) Γ Γ  C    C   δW ext

δW int

where the left and right sides represent the work done by the external force and stress, respectively, δu(x) is a virtual displacement and δe(x) is the virtual Euler–Almansi strain tensor. The displacement vector is defined as the difference between the current and reference configuration u(x) = x(t) − X and the virtual Euler–Almansi strain is given by ) 1( (27) δe(x) = ∇s δu(x) + ∇ sT δu(x) . 2 In curvilinear coordinates, it can be rewritten in covariant form as ) 1( i δeαβ (x) = δu α,β (x) + δu β,α (x) − 2Γαβ (x)δu i (x) (28) 2 i 3 where Γαβ = gi · gα,β is the Christoffel symbol and bαβ = Γαβ is the curvature tensor. Considering the isogeometric concept, the virtual displacement and traction on the membrane are discretized by the T-spline basis functions as δu(x) =

mC ∑

R(x) A δu A

A=1 mC

t(x) =

∑

(29) R(x) A t A

A=1 m

m

C C where m C is the number of global basis functions that defines the capsule and {δu A } A=1 and {t A } A=1 are the virtual mC displacement and traction control vector corresponding to the global basis {R A } A=1 . The left and right terms of Eq. (26) in discrete form are given by ] ∫ mC mC [ ∑ ∑ δW ext = δuTA t B R A (x)R B (x) dΓ (30)

A=1 B=1

Γ

and δW int =

( ) ] ∫ mC [ ∑ 1 1 i δuTA τ αβ (x) R A,β (x)gα + R A,α (x)gβ + R A (x)gα,β − Γαβ R A (x) dΓ 2 2 Γ A=1

(31)

80

J. Maestre et al. / Comput. Methods Appl. Mech. Engrg. 326 (2017) 70–101

respectively. Note that for a given current configuration the Cauchy stress, τ (x), can be evaluated directly using the constitutive equation (Eq. (6)). Eq. (26) can be written in matrix form as St = d

(32) 3m C ×3m C

3m C

where S ∈ ℜ is the global matrix of the system, t ∈ ℜ is the global traction vector on the membrane and d ∈ ℜ3m C is the known right hand side term of Eq. (26). Standard techniques can be used to assemble this system. 3.4. Unperturbed flow We consider a uniform pipe with a fixed cross section (see Fig. 1). Assuming that the axis of the pipe is aligned with the Cartesian x1 -axis, the unperturbed flow through the pipe is given by the Stokes flow equation, ∞ ∞ ∞ µ(v1,22 + v1,33 ) = p,1

(33)

and the continuity equation v1,1 = 0. Note that the transverse components of the velocity field (v2∞ , v3∞ ) are zero as ∞ ∞ well as the pressure gradients ( p,2 , p,3 ) in the section. The flow is constant along the pipe and only depends on the ∞ cross section and longitudinal pressure gradient, p,1 . Analytical solutions are known for some particular section shapes (see for example [116]). For a general shape, we compute the unperturbed flow by the BEM applied to the section of the pipe. To derive a regular BIE, the flow is ∞, p ∞, p divided in two parts: a particular flow, v1 , and a complementary flow, v1∞,c ; such that v1∞ = v1 + v1∞,c . We take the particular flow as, ∞, p

v1

= Ax22 + Bx32 + C x2 x3 + Dx2 + E x3

(34)

∞ that satisfies Eq. (33) for 2(A + B) = p,1 /µ. In this approach, we set the coefficients A = B and the rest of them (C, D, E) are set equal to zero for simplicity. The complementary flow can be computed by the BIE, ∫ ∫ ∂v ∞,c ∂G dl − G 1 dl (35) Cv1∞,c = − v1∞,c ∂n ∂n L ar c L ar c ∞, p

with the boundary condition, v1∞,c = −v1 , on the boundary of the section, L ar c (see Fig. 1). In Eq. (35), G is the bi-dimensional fundamental solution for the Laplace’s problem and C is the jump term [106]. 3.5. Dynamic refinement In dynamic problems with large deformations, the accuracy of the final result is very sensitive to the accumulation of numerical errors in each time step, thus the underlying spatial discretization plays an important role. As the capsule moves and deforms, the mesh resolution must change to capture changes in the solution such as high gradients of the physical variables. To address this issue, we leverage T-spline local refinement. In particular, we focus on geometrically exact h-refinement. Additional details on T-spline h-refinement can be found in [66,101]. To drive the local refinement, we construct a set of error metrics. We measure the size and variation of the tractions, the velocity and the mean curvature (χ = 1/2∇ · n) along the two parametric directions of each T-spline element. The size and curvature metrics allow us to control the geometrical aspects of the model as well as the quality of the mesh through the measurement of the aspect ratio of the T-spline elements. On the other hand, the velocity and traction metrics provide a measure of the quality of the discretization of the physical variables. These metrics are given by ∫ ⎧ ⎪ ⎪ L = dl α ⎪ ⎪ lα ⎪ ⎪ ⏐ ⏐∫ ⏐ ⏐⏐ ⎪ ⎪ ⏐ ⎪ ∆t|lα ⏐ 1 ⏐⏐ ⎪ ⏐ ⎪ ˆ ⎪ ⎨Tα = L α ⏐ dt⏐ = L α l ⏐ ⏐∫ α ⏐ ⏐⏐ (36) ⏐ ⏐ ∆v|lα ⏐ ⎪ 1 ⎪ ⏐ ⏐ ˆ ⎪ Vα = dv = ⎪ ⎪ L α ⏐ lα ⏐ Lα ⎪ ⎪ ⏐ ⎪ ⏐∫ ⏐ ⏐⏐ ⎪ ⏐ ⎪ ⏐ ⏐ ∆χ| 1 ⎪ lα ⎪ ⏐ dχ ⏐ = ⎩χ ˆα = ⏐ ⏐ L L α

lα

α

J. Maestre et al. / Comput. Methods Appl. Mech. Engrg. 326 (2017) 70–101

where lα is the central arclength of the element and α = 1, 2. The elements are candidates for subdivision if ( r e f )κ ( r e f )κ ( r e f )κ } { ˆ ˆ Lα ˆα ref r e f Tα r e f Vα ref χ or < Dr e f L α > min L α , L α , Lα , Lα ˆα ˆα χ ˆα Lβ T V

81

(37)

in which the superscript (·)r e f denotes the set of reference values, Dr e f is the limit aspect ratio and κ is a set order. The adaptive algorithm is based on the management of the aspect ratio and the two central arclengths of the element. The discrepancy between the variation of the physical variables with respect to the reference values acts as a penalty factor reducing the maximum size of the element in areas with high gradients. The element can be divided independently in each parametric direction producing two or four subelements. To avoid excessive refinement, a minimal length (L min α ) is set. Once the elements are selected, they are divided by edge insertion using the T-spline refinement algorithm [101]. It is worthy to mention that the local refinement process is natural in the T-spline context. In contrast to the traditional FEA, the geometry and physical variables are kept unaltered after the refinement since the initial and refined T-spline spaces are nested. However, the local coarsening is not straightforward and the algorithm is complex [117]. To overcome this limitation, new advanced technologies based on Hierarchical T-splines [88] and the B´ezier projection technique [118] have recently appeared. For deformable capsules, two geometries are constructed: one that corresponds to the reference configuration (or unstressed state) and another that corresponds to the current configuration. The evaluation of the mesh quality is performed in the current configuration and the refinement process is carried out in both configurations to ensure the same topology and parametrization. This strategy is general and not all mesh quality measures are required. The selection of the appropriate mesh quality measure depends on the needs of the particular problem. In Section 4 we will demonstrate the effectiveness of this strategy. 3.6. The coupled numerical approach To integrate the system in time, a robust time stepping algorithm is required. In the literature, various schemes have been proposed to simulate deformable capsules in Stokes flows, and most of them are explicit methods [25,28,43,62,98,119,120]. These algorithms are conditionally stable. For our application, an explicit scheme requires unacceptable small time steps to maintain stability and accuracy. Our experience is that instabilities usually manifest when degenerate elements are present in the mesh, when the capsule is very stiff or when the membrane undergoes buckling. To overcome this issue, we use an implicit second-order scheme with an adaptive time step. The Crank–Nicolson scheme can be written as ) ∆tn ( v(tn ) + vk (tn+1 ) k = 1, 2, . . . (38) xk+1 (tn+1 ) = x(tn ) + 2 where tn denotes the current time and k is the sub iteration number. To update the position of the membrane to the next time step, tn+1 , an iterative process is carried out until the normalized L 2 -error of the control points satisfies a √ tolerance error, ∥xk+1 (tn+1 ) − xk (tn+1 )∥ L 2 /(RC m C ) ≤ ϵ, in which RC is the characteristic radius of the capsule (i.e. the radius of the equivalent sphere with the same volume). The time step, ∆tn , is modified in each iteration until the integration scheme converges to the imposed tolerance. In other words, if the number of sub iterations, k, exceeds a maximum number of iteration, Imax , without reaching the tolerance required, or is smaller than a minimum number, Imin , the time step is modify by some factor. For second-order time integration schemes, convergence is expected when [43], Q ∞ ∆tn hCa < O( ) RT π RT3

(39)

where RT is the radius of the tube, h is the representative minimum element size of the capsule (defined as the square root of the element surface) and Ca is the capillarity number. Eq. (39) gives the order of magnitude of the time step needed and it can be used to select the value of the initial time step. In our simulations we set the tolerance error to ϵ = 10−7 . The maximum number of sub iterations is Imax = 8, the minimum number is Imin = 4 and the time scale factors are 0.75 and 1.33, respectively. This selection is made based on empirical criteria and it produces accurate results as shown in the numerical examples (Section 4).

82

J. Maestre et al. / Comput. Methods Appl. Mech. Engrg. 326 (2017) 70–101

The capillarity number and the ratio between the radius of the capsule and the radius of the tube (δ = RC /RT ) are important parameters because they influence the motion and deformation of the capsule. The capillarity number relates the fluid stress on the membrane with the elastic stress. In the case of a capsule flowing through a pipe it is defined as, µvm∞ µQ ∞ = Ca = (40) Gs π RT2 G s where vm∞ is the average velocity flow across the pipe section. Although the BIE satisfies implicitly the null divergence of the velocity field through the fundamental solution, the volume of the capsule can still experience slight changes due to numerical errors. There are three primary sources of numerical errors: (1) errors in the discretization of the surface and velocity field; (2) errors in the numerical integration of the kernels; and (3) errors produced because the time integration scheme is not conservative. These errors are small but can accumulate over time. To control the volume change, we constrain the fluid system to satisfy the condition that the net normal flux across the membrane is zero in each iteration (the interested reader is referred to [62] for more specific details about the modification of the fluid system of equations and the numerical implications of this constriction). The constrained equation is written as, ∫ v(x) · n(x)dΓ = 0 (41) ΓC

For later use, we define the centroid of the capsule at any time by the boundary integral expression: ∫ 1 vc x 2 n i (x)¯xi dΓ x = 2v ΓC i

(42)

where x¯ i is the unit Cartesian vector in the i-direction and v is the actual volume of the capsule. The volume can be written as ∫ x3 n 3 (x)dΓ (43) v= ΓC

Fig. 5 shows the flow chart corresponding to our present numerical procedure. Note that the fluid and mechanical solver are decoupled. Knowing the reference and the current configuration of the membrane at any time step, we apply the PVW (Eq. (26)) to compute the traction on the capsule. Then, introducing the jump stress into the fluid solver system (Eqs. (24), (25) and (41)), we obtain the velocity of the capsule. As pointed out by Pozrikidis [95], the fluid solver is simplified when the internal and external fluids are isoviscous. In this case, the BIE on the capsule and tube decouple and it is only necessary to compute the single layer operator. Once the velocity of the interface is computed, the shape is updated using the integration scheme (Eq. (38)). To avoid constructing an unnecessarily long tube, we relocate the capsule to the center of the tube at each iteration. That is, after updating the capsule, we compute the variation of centroid of the capsule along the longitudinal direction, ∆x1vc = x1vc (tn+1 ) − x1vc (tn ), and subtract this value from the actual position (we assume the tube extends in the x1 -direction). This procedure is repeated until the final time is reached. 4. Numerical examples We solve several numerical examples to show the accuracy and stability of the IGA BEM–FEM approach. As final application of our method, we analyze the behavior of a deformable capsule flowing through a capillary. This problem has important implications in some biomedical flows. For later use, we define the L 2 -error between a known vector field, w, and a discrete solution, wh , as √∫ h ∥w(x) − wh (x)∥2 dΓ (44) ∥w − w ∥ L 2 −Err or = Γ

and the pointwise error as ∥w(x) − wh (x)∥. In the following examples, eight Gauss quadrature points per element along each parametric direction were selected for the numerical evaluation of the integrals. As suggested by Heltai et al. [63], the use of a higher number of quadrature points does not provide a significant increase of the numerical accuracy whereas it penalizes the computational cost.

J. Maestre et al. / Comput. Methods Appl. Mech. Engrg. 326 (2017) 70–101

83

Fig. 5. Flow chart of the numerical procedure. The table summarizes the main equations used during the procedure.

Fig. 6. Rigid Spheroid moving in a quiescent fluid.

4.1. Uniform flow past a spheroid Our first numerical example is a constant flow over a fixed spheroid. The aim of this example is to analyze the spatial convergence of the method. The geometry is similar to a deformed capsule and we can study the impact of the geometry on the accuracy of the method. This problem has an analytical solution for the hydrodynamic forces on the spheroid (the analytical expressions can be found in [121]). The reduced BIE for this problem is: ˆ Γ v)(x0 ) − (G Γ ∆f)(x0 ) + v∞ (x0 ) v(x0 ) = (H

(45)

We can assume that either the fluid is moving with respect to a fixed spheroid (v = 0, v∞ ̸= 0) or that the sphere is moving with respect to the fluid (v ̸= 0, v∞ = 0). This change of observer leads to two different numerical formulations. We assume the second option and check the accuracy of the single and double layer kernels. We define the spheroid through its volume V = 4/3πa 2 b = 1 and the relationship between the semi-major axis and the semi-minor axis, b/a = 1.5 (see Fig. 6), such that the semi-major axis is oriented vertically (x3 -direction). The incidental angle of the flow is set to θ = π/4 with respect to the horizontal plane (x1 x2 -plane), thus the velocity of the body is given by v = −(cos(θ )¯x1 + sin(θ )¯x3 ) where x¯ i are the Cartesian vectors. We study two different mesh constructions of the spheroid (Fig. 7). The first is a structured mesh with the poles on the semi-major axis (Fig. 7(a)). The second (Fig. 7(b)) is an unstructured mesh without poles. The model is constructed by transforming a regular cube to a spheroid. In both cases, we set the weights to be one and we use smooth basis functions (expect at the poles where the continuity is C 0 ). Four different discretizations of each model are used. For

84

J. Maestre et al. / Comput. Methods Appl. Mech. Engrg. 326 (2017) 70–101

Fig. 7. (a) Structured mesh and (b) unstructured mesh of a spheroid. In (a) the red spherical control points denote the poles and in (b) the red triangular points are extraordinary points.

Fig. 8. L 2 -Error of the discretized normal field for the structured and unstructured mesh.

the structured mesh, we set 42, 146, 546 and 2114 control points and for the unstructured mesh 56, 218, 866 and 3458 control points. The coarsest mesh is shown in Fig. 7. During the refinement each B´ezier element is subdivided into four elements. To achieve accurate solutions, an accurate discretization of the normal field is required. The traction field depends on the normal field through the relation f = τ n. The analytical normal field cannot normally be represented exactly by the T-spline basis functions, introducing errors into the method. Consequently, the error of the discretized traction vector, without considering other numerical aspects, is governed by the error of the discretized normal vector. In Fig. 8, we show the L 2 -error between the analytical normal field and the discretized normal field (computed using Eq. (2)) for the structured and unstructured mesh as a function of the number of control points. The errors converge quasi-linearly when a logarithmic scale is used. The average convergence ratio, defined as the logarithmic slope, is 2.03 for the structured mesh and 1.62 for the unstructured mesh. Note that the structured mesh contains degenerate elements. These degeneracies reduce the accuracy of the method if the degenerate elements are not treated. Taus et al. [102] propose a reparametrization based on the singular-value decomposition theorem accompanied by an increase in the number of quadrature points to improve the numerical integration of the kernels. On the other hand, Heltai et al. [63] implemented a non-singular BIE providing optimal results. This approach is based on removing the singularity of the single and double layer operators, obtaining a well conditioned system of equation. This technique requires a smooth geometry and assumes that the normal is unique

J. Maestre et al. / Comput. Methods Appl. Mech. Engrg. 326 (2017) 70–101

85

Fig. 9. (a) L 2 -error of the pressure and (b) percentage error of the drag force on the spheroid for the unstructured and the structured mesh using different techniques. Table 1 Errors of the normal vector, pressure and drag force on the spheroid for the structured mesh, discretized by 546 control points, and using the Bézier transformation technique. For each subdivision level, the original degenerate Bézier element is subdivide in 2n Bézier elements in the azimuthal direction. Subdivision level, n Normal vector (L 2 -error) Pressure (L 2 -error) Drag (%Error)

0 1.332 · 10−3 6.499 · 10−4 1.735 · 10−6

1 3.494 · 10−4 2.362 · 10−4 1.682 · 10−6

2 8.947 · 10−5 1.825 · 10−4 1.679 · 10−6

3 8.947 · 10−5 1.781 · 10−4 1.679 · 10−6

at the collocation points. In the presence of degenerate elements, the collocation points are often placed at the poles where the continuity of the surface is locally C 0 . For a fixed spheroid the normal is unique at these points though it cannot be directly computed from the underlying parametrization. However, in deformable capsules the normal can become non-unique and the method is not applicable. Fig. 9 shows the results obtained with the standard BIE for the structured and the unstructured mesh, and those obtained by the desingularization of the kernels and by the B´ezier transformation technique (explained in Section 3.1) for the structured mesh. For the B´ezier transformation, the original degenerate elements are subdivided into eight B´ezier elements in the azimuthal direction to provide an accurate approximation of the original B´ezier elements. Fig. 9(a) depicts the L 2 -error of the normal pressure on the spheroid. The error obtained with the standard formulation for the structured mesh is almost constant for all discretizations. This is because the L 2 -error is dominated by the error in the degenerate elements. For the remaining cases, the error rapidly converges regardless of whether the desingularization technique or the B´ezier transformation are used. The average convergence ratios are 1.61 for the structured mesh (with the B´ezier transformation technique) and 1.64 for the unstructured mesh. The total drag error on the spheroid is shown in Fig. 9(b). Erratic behavior is observed for the standard formulation and the structured mesh. The other cases converge almost linearly (in a logarithm scale). The rates are 2.21 for the structured mesh and 2.28 for the unstructured mesh. In general, the structured mesh exhibits smaller errors than the unstructured mesh when the singular elements are treated correctly, but the converge rates are similar and they are on the same order as the converge rates of the normal vectors. In Table 1 we summarize the errors of the normal vector, pressure and drag on the spheroid for the structured mesh (discretized with 546 control points), using the B´ezier transformation and various subdivision levels. These errors are reduced with the subdivision of the original element until reaching an almost constant value from the level 2. Similar behaviors are found for other mesh discretizations. Fig. 10 shows the point-wise error of the pressure on the spheroid surface. The error is concentrated around the poles when no reparametrization is performed. The error is smaller and more uniformly distributed after the

86

J. Maestre et al. / Comput. Methods Appl. Mech. Engrg. 326 (2017) 70–101

Fig. 10. Point-wise error of the pressure on the spheroid for (a) the structured mesh (546 control points) without the reparametrization, (b) the structured mesh with the Bézier transformation technique and (c) the unstructured (866 control points).

application of the B´ezier transformation technique to the degenerate elements. For the unstructured mesh, the error is also distributed throughout the surface of the spheroid, but it is higher around the extraordinary points probably because the constrained approach applied to force the continuity of the basis functions across these points. 4.2. Capsule in shear flow We now analyze the temporal stability and accuracy of the coupled IGA BEM–FEM formulation by simulating the deformation of a capsule immersed in a shear flow. Although this problem has no known analytical solution, it has been studied extensively in the literature [25,43]. We focus on the case of an initial spherical isoviscous (λ = 1) capsule with a Neo-Hookean law for the membrane. The BIE for this problem is v(x0 ) = −(G Γ ∆f)(x0 ) + v∞ (x0 )

(46)

where the unperturbed flow is given by v = γ˙ x3 x¯ 1 and γ˙ is the shear rate. The flow deforms the capsule until reaching a steady shape. At this constant shape, the membrane moves and rotates with respect to the internal fluid. This phenomenon is known as tank-treading. In this state, the capsule shape is similar to an ellipsoid with the major-axis contained in the shear plane. Consequently, the deformation of the capsule is usually quantified through the Taylor parameter, ∞

L1 − L2 (47) L1 + L2 where L 1 and L 2 are the longest and shortest semi-axes measured from the centroid to the surface, respectively. For this problem, this parameter provides a measure of the change of shape. The motion and deformation of the capsule is characterized by the capillarity number which, for this case, is defined as Ca = µγ˙ RC /G s . There are low and high critical capillarity numbers for which the capsule will not reach a stable equilibrium. We compare the temporal deformation of the capsule obtained by our formulation with those presented by Lac et al. [25] for various stable values of capillarity, Ca = 0.3, 0.45 and 0.6, and different meshes. Without loss of generality, we consider an initial sphere with unit radius. As in the previous example, we consider two topological representations of the capsule. The first is a regular structured mesh and the other is an unstructured mesh. The structured mesh has 842 control points, and the unstructured mesh, constructed through a cube projection, has 866 control points. The B´ezier transformation technique is used to reparameterize the degenerate elements. The initial dimensionless time step is set to γ˙ ∆t0 = 0.01Ca. The simulation is stopped at the dimensionless time γ˙ tend = 25Ca to guarantee that the steady state is reached. In Fig. 11(a) we show the time evolution of the Taylor parameter for the unstructured mesh and three different capillarity numbers. The computed results are very close to those reported in the literature [25] for all capillarity numbers considered. The angles of the principal direction of the capsule with respect to the horizontal at the steady state, θ , (Fig. 11(b)) are also measured and are in agreement with those obtained in [25]. Similar results are found for the structured mesh. The relative discrepancy of the deformation parameter and orientation of the capsule between the unstructured and the structured mesh are practically negligible, lower than O(10−3 )% and O(10−2 )%, respectively. D12 =

87

J. Maestre et al. / Comput. Methods Appl. Mech. Engrg. 326 (2017) 70–101

Fig. 11. (a) Evolution of the Taylor deformation parameter of the capsule (Neo-Hookean membrane) in a shear flow for the unstructured mesh and three capillarity numbers (Ca = 0.3, 0.45 and 0.6). (b) Steady shape of the capsule in the shear plane for the unstructured mesh. The inserted table contains the angles of the principal axis of the deformed capsule. Table 2 Average and maximum volumetric error of the capsule in shear flow during the simulation for both the structured and unstructured mesh and various capillarity numbers (Ca = 0.3, 0.45 and 0.6). Ev =

∆v V %

Ca = 0.3

Unstruc. Struc.

Volume constrain Non volume constrain Volume constrain Non volume constrain

Ca = 0.45

Ca = 0.6

Avg(E V )

max(E V )

Avg(E V )

Max(E V )

Avg(E V )

Max(E V )

1.922 · 10−4 4.449 · 10−4 3.326 · 10−5 2.815 · 10−4

4.285 · 10−4 9.892 · 10−4 4.580 · 10−5 5.794 · 10−4

1.443 · 10−4 1.443 · 10−3 6.470 · 10−5 6.582 · 10−4

7.603 · 10−4 3.358 · 10−3 1.031 · 10−4 2.717 · 10−3

2.065 · 10−4 2.957E · 10−3 7.588 · 10−5 1.671 · 10−3

1.173 · 10−3 5.841E · 10−3 1.547 · 10−4 6.324 · 10−3

Table 3 Average time step of the capsule in a shear flow during the simulation for both the structured and unstructured mesh and various capillarity numbers (Ca = 0.3, 0.45 and 0.6). Avg(∆t · γ˙ /Ca) Unstructured Structured

Ca = 0.3

Ca = 0.45

Ca = 0.6

4.480 · 10−2

4.299 · 10−2

1.033 · 10−2

1.100 · 10−2

4.094 · 10−2 1.16510−2

An illustrative indicator to verify the accuracy of the method is the time evolution of the volumetric change of the capsule (E v (t) = ∆v/V). This indicator captures the effect of the volumetric constrain (Eq. (41)). The results are summarized in Table 2. It is observed that the volumetric error remains low for all capillarities and both the structured and unstructured mesh (Avg(E v ) ≤ O(10−3 )%). However, the volumetric error for the structured mesh is smaller because it requires a smaller time step compared with the unstructured mesh (see Table 3). Moreover, this error increases with the capillarity, as it is expected because the capsule suffers larger deformations, and significantly decreases when the volumetric constrain equation is used. Note that the structured mesh contains elements whose aspect ratio (ratio between length and width) is highly variable. In fact, the elements around the poles show an aspect ratio far from the ideal value of 1. This high aspect ratio coupled with the loss of smoothness and the singularity of the degenerate elements, produces an important source of numerical errors. Joneidi et al. [62] found numerical instabilities around the poles for a deformable capsule discretized with a NURBS. Similar finding were reported by Lac et al. [25] using bi-cubic B-spline functions. In [25] the authors analyzed the influence of the initial position of the poles and found that the simulations were unstable

88

J. Maestre et al. / Comput. Methods Appl. Mech. Engrg. 326 (2017) 70–101

Fig. 12. (a) Instability of the capsule around the poles located on the x2 -axis when the membrane equilibrium equation in strong form (Eq. (3)) is used (at γ˙ t = 0.664), and (b) steady state of the capsule using the variational equation (Eq. (26)) for the Capillarity number Ca = 0.45. The transparent sphere corresponds to the initial state and the solid surface is the deformed state.

when they were placed on axes which are out of the shear plane (x1 x3 -plane) where the membrane is compressed. The results showed (for the structured mesh) was obtained for the poles contained on the shear plane. However, we find that these instabilities are avoided by reparameterizing the degenerate elements, applying the B´ezier transformation technique (or the technique proposed by Taus et al. [102]) and using an adaptive implicit temporal integration scheme. Additionally, we carried out simulations with the poles contained in the x2 -axis and for an intermediate capillarity number, Ca = 0.45. Using a fixed time step (∆tn = ∆t0 = 0.01Ca), the temporal integration scheme does not converge at the early stages and it must be reduced by at least 25% to achieve stable results. This is overcome using the proposed adaptive time integration scheme. Another important aspect affecting the overall stability of the method is how the membrane load is computed. Our simulations reveal that instabilities appear when the equilibrium equation in strong form (Eq. (3)) is employed, even when an adaptive time integration scheme is used. The variational equation (Eq. (26)) provides stable results independently of the position of the poles. Fig. 12 shows a typical instability produced around the poles. 4.3. Dynamic mesh In this section we show the impact of dynamic mesh refinement on the solution. Two examples are shown. The first example is the movement of a micro-drop immersed in a four-roll flow. Although it is not a capsule, this example is significant because the error is very sensitive to the quality of the mesh. We analyze a particular case for which numerical and experimental results have been previously reported in the literature [120,122]. Following the notation from [120], the unperturbed flow can be expressed as: ⎧ ) ( ) ) γ˙ (( ⎪ ⎪ ⎪ ⎨v1 = 2 1 + a f x1 + 1 − a f x3 v2 = 0 (48) ⎪ ) ( ) ) ⎪ γ˙ (( ⎪ ⎩v3 = −1 + a f x1 − 1 + a f x3 2 where γ˙ is the shear rate and a f is a parameter which modify the flow direction and it is set to 0.6. The viscosity ratio C is λ = 0.118 and the capillarity number is given by Ca = γ˙ µR = 0.196, where σ the superficial tension and RC the σ radius of the drop (see Fig. 13). The computation of the jump stress on the drop surface depends only on geometric aspects (i.e. ∆f = (∇ · n)n). Consequently, the motion of the interface is described by: ˆ Γ v)(x0 ) − (G Γ ∆f)(x0 ) + v∞ (x0 ) v(x0 ) = (1 − λ)(H

(49)

By the action of flow, the drop is deformed into a shape similar to an ellipsoid in the steady regime. Note that the fluid particles on the interface are not bounded. In other words, the fluid particles can move freely over the interface due to the net fluid stress on both sides. This can produce a migration of the fluid particles to certain zones of the

J. Maestre et al. / Comput. Methods Appl. Mech. Engrg. 326 (2017) 70–101

89

Fig. 13. Unperturbed velocity field for a four-roll flow with a f = 0.6.

interface. Coupled with the large deformations of the drop, this can lead to large distortions of the mesh. To manage mesh distortions we use dynamic local refinement. The simulation starts from a sphere with unit radius. A structured mesh with 222 control points is used. This initial discretization is not fine enough to capture the large deformation of the drop. Given the characteristics of this problem, the suitability of the mesh is evaluated as a function of the variation of the curvature and the aspect ratio of the elements. It should be noted that the tractions are proportional to the mean curvature, thus this parameter is representative of the geometry discretization as well as the discretization of the physical variables (traction and ref ref velocity). The reference parameters are set to: L α = 0.32RC , χ ˆα = 2/RC2 , Dr e f = 2, κ = 2 and L min = 0.05RC α where α = 1, 2. In the degenerate elements, refinement is not applied. The initial time step is γ˙ ∆t = 0.005. The simulation ends once the steady state is reached at γ˙ tend = 10. Fig. 14 shows the time evolution of the Taylor parameter. The present results are very close to the numerical results ∞ reported in Wang and Dimitrakopoulos [120]. Moreover, the deformation of the drop at the steady state (D12 = 0.381) ∞ is in good agreement with the experimental findings of Bentley and Leal [122] (D12 = 0.381). Additionally, we simulate the deformation of a drop with a fixed mesh. If a refinement is not applied, the accuracy of the results deteriorates significantly (the relative discrepancy is 8.312%). This reveals the importance of a suitable discretization of the drop. In Fig. 15 the discretization of the drop is shown at three different stages. As the drop deforms, the variation of the curvature increases rapidly around the two endpoints. This corresponds to the principal inertial axis of the drop. To capture this variation, various local refinements are performed. At the end of the simulation the drop contains a total of 844 control points. We then studied the movement of an isoviscous capsule in a hyperbolic flow following a Skalak law with C = 1. The unperturbed flow is expressed by Eq. (48) with a f = 1. The capsule again evolves into a shape similar to an ellipsoid. This problem has been analyzed by several authors using a BEM approach [25,28,43]. Lac et al. [25] found numerical instabilities when the capillarity number is larger that the value 0.7. Dodson and Dimitrakopoulus [28] and Walter et al. [43] computed stable results for higher values of capillarity number (Ca > 0.7). They observed that the capsule acquires very high mean curvature during the transitory phase which can lead to numerical problems. Given these findings, we carry out some simulations for capillarity numbers Ca = 0.45, 0.60 and 0.8. The use of a dynamic local refinement scheme is required, otherwise an excessively fine discretization is needed during the entire simulation. The simulation is initialized using sphere with unit radius. The sphere is discretized with an unstructured mesh with 218 control points (Fig. 16(a)). Mesh distortion and refinement are managed through the variation of the traction ref ˆαr e f = 8G s /RC Dr e f = 2.5, L min on the membrane. The reference values are set to: L α = 0.35RC , T = 0.035RC α and κ = 1 where α = 1, 2. These parameters were established experimentally and are based on the deformation of the capsule at the steady state. We find stable results for all capillarities and Taylor parameters are in very good agreement with those reported by Dodson and Dimitrakopoulos [28] and Walter et al. [43] (see Table 4) at the steady state. Dodson and Dimitrakopoulos [28] employed a spectral BEM model with a range of 216–1176 spectral points and used a smooth

90

J. Maestre et al. / Comput. Methods Appl. Mech. Engrg. 326 (2017) 70–101

Fig. 14. Time evolution of the Taylor parameter for a drop in a four-roll flow for Ca = 0.196 λ = 0.118. The circles denote the times at which a refinement is applied to the drop. The dash marker ‘-’ on the right corresponds to the Taylor deformation experimentally found by Bentley and Leal [122] at the steady state. The inset (a) is the experimental shape reported in [122] (reproduced with permission) and (b) shows the shape of the drop obtained by the present approach at the steady state. Table 4 Taylor parameter of a isoviscous Skalak capsule (C = 1) in a hyperbolic flow at the steady state for different Capillarities. ∞ Taylor parameter D12

Present results Dodson and Dimitrakopoulos [28] Walter et al. [43]

Ca = 0.45

Ca = 0.6

Ca = 0.8

0.548 0.548 0.55

0.595 0.595 0.60

0.640 0.64

technique between spectral elements, and Walter et al. [43] used a standard BEM–FEM model based on quadratic elements with 2562 points for the computations. The discretization of the capsule in the steady state as well as the variation of the tractions are illustrated in Fig. 16. The capsule suffers very large deformations (larger for higher values of Ca) and the variation of tractions is very high around the tips of the capsule. This requires a fine mesh in these areas. The total number of control points at the steady state is 706 for Ca = 0.45, 906 for Ca = 0.6 and 1478 for Ca = 0.8. The volumetric error remains small, O(10−3 %), for all simulations. We emphasize that our approach is stable and accurate without applying smoothing techniques. Fig. 17 shows a comparison of the capsule shape for Ca = 0.80 at the steady state between a fixed mesh and a locally refined mesh. Despite the fact that the fixed mesh has a low number of elements, the capsule shape for both meshes is surprising similar. The maximum discrepancy occurs at the tips due to the high variations in the tractions and curvature which cannot be accurately captured by the coarse mesh. The error in the tractions between the fixed and locally refined mesh is high (17.32%) as shown in Fig. 18. 4.4. Capsule in a capillary We now analyze a biological flow consisting of a deformable capsule moving in a capillary. The capillary is a long tube with a constant cross section. The section is taken from a real picture corresponding to a brain capillary section of a rat (see Fig. 11b in [123]). The capsule is a spherical Skalak membrane (with C = 1) that contains an isoviscous fluid. The capsule moves freely and is deformed by the action of a constant flow across the capillary which perturbs the behavior of the capillary flow.

J. Maestre et al. / Comput. Methods Appl. Mech. Engrg. 326 (2017) 70–101

91

Fig. 15. Discretization of the drop in a four-roller flow applying a dynamic refinement (a) at the beginning γ˙ t = 0, (b) in an intermediate state γ˙ t = 0.509 and (c) in the steady state γ˙ t = 10. (d), (e) and (f) variation of the mean curvature of the vertical section of the drop (x1 x3 -plane) represented in a polar coordinate system at the three respective times. Dashed line: scale representation of the section of the drop (E = 3 : 1); solid line: variation of the curvature ( dχ χ ). dl ); and dotted line: reference value of the variation of the curvature (ˆ

Fig. 19 shows an illustration of the model used for the simulation. For convenience, the two main inertial axes of the capillary section are aligned with the Cartesian x2 and x3 axes such that centerline of the tube, defined as the centroid of the cross section, coincides with the x1 -axis. The tube has a length of L T = 2π RT which is long enough as it was assumed in Section 3.2 [95]. The ratio between the capsule radius and the equivalent tube radius is set to δ = RC /RT = 0.6. The center of the capsule is placed in the middle of the tube and at a radial distance rC = 0.2RT under the centerline. The capsule is discretized using an unstructured mesh with 602 control points and the tube is modeled with a structured mesh with 21 × 20 control points. Three levels of refinement are applied along the tube, concentrating the mesh around the capsule. We considered representative values of the shear modulus of a biological capsule, G s = 6 · 10−6 N/m, the fluid viscosity, µ = 1.2 · 10−3 Ns/m2 , and the average unperturbed flow velocity, vm∞ = 0.5 mm/s. In this case, the capillarity number is equal to Ca = 0.1 (Eq. (40)). In Fig. 20 we show the unperturbed flow velocity across the capillary section. For comparative purposes, we also simulate the capsule flowing in a cylindrical tube whose flow is given by the Poiseuille law. Additionally, we assume a positive pressure jump ( pos ) between the interior and exterior of the capsule. This situation occurs naturally in capsules with semi-permeable membranes because of different concentrations of any substance on either side of the membrane (i.e., osmotic pressure). For a spherical isotropic capsule, this pressure jump pi pi can be taken into account through a preinflation of the capsule from a radius RC to a radius of RC = (1 + α)RC , causing an isotropic internal stress τiso = 1/2 pos RC [24–26]. The capsule preinflation helps to prevent the buckling of the membrane. Buckling phenomena cannot be handled by the present approach due to the lack of bending stiffness in the membrane. We assume a low expansion coefficient, α = 0.02.

92

J. Maestre et al. / Comput. Methods Appl. Mech. Engrg. 326 (2017) 70–101

Fig. 16. Locally refined capsule in a hyperbolic flow (a) at the initial time and at the steady state for (b) Ca = 0.45, (c) C = 0.60 and (d) Ca = 0.8. (e), (f) and (g) show the variation in the tractions in the meridional plane of the capsule (x1 x3 -plane) represented in a polar coordinate system at the three respective capillarity numbers. Dashed line: scale representation of the section of the drop (E = 30 : 1); solid line: variation of the tractions ˆ ( dt dl ); and dotted line: reference value of the variation of the traction ( T ).

Fig. 17. Capsule shape in a hyperbolic flow with Ca = 0.8 at the steady state. Solid line: locally refined mesh; dashed line: fixed mesh.

Fig. 21(a) shows the time evolution of the capsule shape flowing through the tube. The capsule is deformed rapidly and moves in the direction of the flow as well as it ascend vertically (this phenomenon is known as migration). In Fig. 21(b) we show a longitudinal section of the capsule (x1 x3 -plane) for different time steps. The frame of reference

J. Maestre et al. / Comput. Methods Appl. Mech. Engrg. 326 (2017) 70–101

93

Fig. 18. Tractions on the membrane for a capsule in a hyperbolic flow for Ca = 0.8 at the steady state. (a) Locally refined solution and (b) fixed mesh solution.

Fig. 19. T-spline model for a capsule flowing in a capillary: (a) transversal section view and (b) perspective view of the model. For clarity, a vertical cut of the perspective view is shown. The main inertial axes of the capillary section are coincident with the Cartesian axes and the centroid is placed in the origin.

∞. Fig. 20. Unperturbed flow velocity across the capillary, v1∞ (x)/vm

is moved with the centroid of the capsule along the flow direction. As shown, the capsule is rapidly deformed and the back part of the capsule acquires a convex shape. As the capsule ascends, it tends to have a parachute shape and the transversal section of the capsule (x2 x3 -plane) becomes quasi-elliptical, similar to the tube shape (see Fig. 21(c)).

94

J. Maestre et al. / Comput. Methods Appl. Mech. Engrg. 326 (2017) 70–101

Fig. 21. (a) Three-dimensional representation of the capsule shape flowing in a tube for different time steps. (b) Longitudinal section (x1 x3 -plane) and (c) transversal section (x2 x3 -plane) of the capsule for different time steps, respectively. The white lines around the capsule, corresponding to the initial three-dimensional view of the capsule, show the section marks of the capsule.

We also analyzed the motion of the capsule along the tube. In Fig. 22(a) and b we show the trajectory of the centroid of the capsule projected on the transversal plane and the transversal velocity (or migration velocity), vC M for all time, respectively. It is observed that the capsule asymptotically migrates toward a point of maximum flow velocity, reducing the migration velocity with the distance to this point. The maximum migration velocity of the capsule in the capillary is less pronounced than for the cylindrical tube and the migration velocity decays more slowly with time. Note that this phenomenon does not happen in the case of a neutral buoyant solid sphere where the particle moves parallel to the tube-axis. The migration phenomenon is due to the deformability of the capsule [29,95,124–129]. Fig. 22(c) shows the time evolution of the horizontal velocity of the capsule, vC H . After a transient phase, the capsule tends to a constant velocity as it goes toward the point of maximum velocity. This velocity is slightly lower for the capillary tube than for cylindrical tube. Moreover, a tank-treading effect is found because of the torque exercised by the flow on the capsule (see Fig. 22(d)). The membrane has a net tangential velocity (VC T ) that produces a clockwise rotation. This effect is reduced with time, decaying more slowly for the capillary than for the cylindrical tube due to a slower migration phenomenon. As is well known, the presence of the capsule in the tube causes an increase of the flow resistance that results in an additional pressure drop (known as disturbed pressure drop), ∆ p D . In Fig. 23 we show the temporal disturbed pressure drop in the tube. The disturbed pressure drop increases rapidly from zero until a maximum peak and then tends to a constant value as the capsule migrates to the point of maximum velocity. The capillary offers a higher resistance than the cylindrical tube because of the section shape and capsule deformation. In fact, the unperturbed pressure drop in the capillary is also higher. The additional pressure drop is an important parameter in biological flows since it is related to the apparent viscosity of the flow.

J. Maestre et al. / Comput. Methods Appl. Mech. Engrg. 326 (2017) 70–101

95

Fig. 22. (a) Trajectory of the centroid of the capsule projected in the transversal plane (x2 x3 -plane). The isolines denotes the flow across the capillary and the cross marker denote the point of maximum velocity. (b) Time evolution of the migration velocity of a capsule on the transversal plane. (c) Time evolution of the horizontal velocity of a capsule in the flow direction.(d) Time evolution of the average tangential velocity in the ∞ /R = 5.00. meridian plane (x1 x3 -plane) of the capsule. The inset (e) shows the tangential velocity distribution of the capsule at the time tvm T

Fig. 23. Time evolution of the disturbed pressure drop due to the presence of a capsule flowing in the tube.

5. Conclusions In this paper, we have presented a three dimensional isogeometric coupled BE–FE analysis in the time domain for the study of deformable capsules. The approach has been developed using analysis suitable T-spline basis functions. The T-spline basis provides a high order of approximation and continuity for both the geometry and physical variables. This continuity and smoothness of the basis are an important property since it confers stability and robustness on the method. Moreover, the high order of the basis provides the high differentiability needed for the accurate computation of the jump stress.

96

J. Maestre et al. / Comput. Methods Appl. Mech. Engrg. 326 (2017) 70–101

In the numerical examples, we examined the spatial accuracy of the present approach analyzing a constant flow past a spheroid. We studied two different topological discretizations: a structured mesh and an unstructured mesh. The structured mesh contains degenerate elements that are a source of numerical problems. We showed that by applying a suitable treatment of degenerate elements both meshes properly converge with nearly the same order of accuracy. We also evaluated the temporal stability and accuracy of the method by simulating a deformable capsule in a shear flow. It was shown that numerical instabilities can appear in the capsule when degenerate elements are present and the membrane equilibrium is applied in strong form. The use of a weak formulation, based on the principle of virtual work, helps to stabilize the method and it reduces the continuity requirement of the basis functions. This linked with the application of an adaptive implicit temporal integration algorithm leads a more robust approach. We compared our results with the literature and we found a good agreement for both structured and unstructured meshes. However, this last mesh is advantageous since it does not contain degenerate elements and the aspect ratio of the elements is more homogeneous. In addition, exploiting the local h-refinement property of the T-spline, we showed the ability of the method to dynamically adapt the mesh of the capsule to the need of the specific problem. This technique is especially important for cases in which the membrane suffers large deformations. Finally, we applied the method to the simulation of a deformable capsule in a capillary to show its potential use in complex biological flows. Acknowledgments This study has been supported by the Spanish Ministerio de Econom´ıa y Competitividad under projects CTQ201346799-C2-1-P and DPI2016-75791-C2-1-P. M. A. Scott was supported by AFOSR Grant FA9550-14-1-0113. These supports are gratefully acknowledged. Appendix. Transformation from degenerate Bézier patch to non-degenerate In Fig. A.24 a schematic illustration of the transformation process from the degenerate quadrilateral B´ezier patch to the non-degenerate B´ezier patches is shown. Using the homogeneous coordinate representation, the transformation from the original degenerate B´ezier patch to the triangular B´ezier patch is given by, ⎤ Qw ⎢ 00 ⎥ ⎢Qw ⎥ ⎢ 10 ⎥ ⎢ w⎥ ⎢Q20 ⎥ ⎢ w⎥ ⎢Q ⎥ ⎥ ⎤ ⎢ 30 ⎥ 0 ⎢ ⎥ ⎢Qw 01 ⎢ 0 ⎥ ⎢Qw ⎥ ⎥ ⎢ 11 ⎥ 0 ⎥⎢ w ⎥ ⎥ Q ⎥ 1/4⎥ ⎢ 21 ⎥ ⎥ ⎥⎢ ⎥ ⎢ 0 ⎥ ⎢Qw ⎥ ⎥ ⎢ 31 0 ⎥ ⎢Qw ⎥ 02 ⎥ ⎥ 0 ⎥ ⎥⎢ ⎢Qw ⎥ ⎥ 12 ⎥ 0 ⎥⎢ ⎢ w⎥ Q ⎥ 0 ⎦⎢ ⎢ 22 ⎥ 0 ⎢Qw ⎥ ⎢ 32 ⎥ ⎢ w⎥ ⎢Q03 ⎥ ⎢ w⎥ ⎢Q ⎥ ⎢ 13 ⎥ ⎢ w⎥ ⎣Q23 ⎦ Qw 33 ⎡

⎤ ⎡ Tw 003 ⎢ w ⎥ ⎡ ⎢T ⎥ 1 ⎢ 102 ⎥ ⎢ w ⎥ ⎢0 ⎢T201 ⎥ ⎢ ⎢ w ⎥ ⎢0 ⎢T ⎥ ⎢ ⎢ 300 ⎥ ⎢0 ⎢ w ⎥ ⎢ ⎢T012 ⎥ ⎢0 ⎢ w ⎥=⎢ ⎢T ⎥ ⎢0 ⎢ 111 ⎥ ⎢0 ⎢ w ⎥ ⎢ ⎢T210 ⎥ ⎢0 ⎢ w ⎥ ⎢ ⎢T ⎥ ⎣0 ⎢ 021 ⎥ ⎢ w ⎥ 0 ⎣T120 ⎦ Tw 030

0 1 0 0 0 −1/4 0 0 0 0

0 0 1 0 0 0 0 0 0 0

0 0 0 1/4 0 0 0 0 0 0

0 0 0 0 1 0 0 0 0 0

0 0 0 0 0 3/4 0 0 0 0

0 0 0 0 0 0 0 0 0 0

0 0 0 1/4 0 0 0 0 0 0

0 0 0 0 0 0 0 1 0 0

0 0 0 0 0 3/4 0 0 0 0

0 0 0 0 0 0 0 0 0 0

0 0 0 1/4 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 1

0 0 0 0 0 −1/4 0 0 1 0

0 0 0 0 0 0 1 0 0 0

(A.1)

In Eq. (A.1) we have been assumed the degeneration is located in the right side of the quadrilateral B´ezier patch. The transformation for other localization of the degenerate side is straightforward by applying a rotation to the conversion matrix. The subdivision of the triangular B´ezier patch in three quadrilateral patches depends on the choice of the points {D0 , D1 , D2 , D3 } that define the domain of each quadrilateral patch (see Fig. A.24). In this work we locate the points {D1 , D2 , D3 } in the center of each side of the triangular B´ezier domain, and the internal point D0 in the barycenter.

J. Maestre et al. / Comput. Methods Appl. Mech. Engrg. 326 (2017) 70–101

97

Fig. A.24. (a) Degenerate cubic quadrilateral Bézier patch, (b) cubic triangular Bézier patch and (c) set of non-degenerate cubic quadrilateral Bézier patches.

This usual configuration provides a well-balanced partition of the triangular B´ezier patch. The computation of each quadrilateral partition is as follows [130]: ⎤ ⎡ ˆw Q 1,00 ⎥ ⎢ w ⎥ ⎢Q ⎢ˆ1,10 ⎥ ⎢ w ⎥ ⎡ ˆ ⎥ ⎢Q 1 ⎢ 1,20 ⎥ ⎢ˆw ⎥ ⎢ 1/2 ⎢Q1,30 ⎥ ⎢ ⎥ ⎢ 1/4 ⎢ w ⎥ ⎢ ⎢Q ⎢ˆ1,01 ⎥ ⎢ 1/8 ⎢ w ⎥ ⎢ ˆ ⎥ ⎢ 1/2 ⎢Q ⎢ 1,11 ⎥ ⎢ ⎢ˆw ⎥ ⎢ 5/18 ⎢Q1,21 ⎥ ⎢ ⎥ ⎢11/72 ⎢ ⎢Q ⎢ w ⎥ ⎢ˆ1,31 ⎥ ⎢ 1/12 ⎢ w ⎥=⎢ ⎥ ⎢ 1/4 ⎢Q ˆ ⎢ 1,02 ⎥ ⎢11/72 ⎢ˆw ⎥ ⎢ ⎢Q1,12 ⎥ ⎢ 5/54 ⎥ ⎢ ⎢ ⎢ ⎢Q w ⎥ ⎢ˆ1,22 ⎥ ⎢ 1/18 ⎢ w ⎥ ⎢ 1/8 ⎥ ⎢ ⎢Q ˆ ⎢ 1,32 ⎥ ⎢ 1/12 ⎢ˆw ⎥ ⎢ ⎢Q1,03 ⎥ ⎣ 1/18 ⎥ ⎢ ⎢ˆw ⎥ 1/27 ⎢Q1,13 ⎥ ⎥ ⎢ w ⎥ ⎢Q ˆ ⎣ 1,23 ⎦ ˆw Q 1,33

0 1/2 1/2 3/8 0 5/18 22/72 3/12 0 11/72 10/54 3/18 0 1/12 2/18 3/27

0 0 1/4 3/8 0 0 11/72 3/12 0 0 5/54 3/18 0 0 1/18 3/27

0 0 0 1/8 0 0 0 1/12 0 0 0 1/18 0 0 0 1/27

0 0 0 0 1/2 5/18 11/72 1/12 2/4 22/72 10/54 2/18 3/8 3/12 3/18 3/27

0 0 0 0 0 3/18 14/72 2/12 0 14/72 13/54 4/18 0 2/12 4/18 6/27

0 0 0 0 0 0 3/72 1/12 0 0 3/54 2/18 0 0 1/18 3/27

0 0 0 0 0 0 0 0 1/4 11/72 5/54 1/18 3/8 3/12 3/18 3/27

0 0 0 0 0 0 0 0 0 3/72 3/54 1/18 0 1/12 2/18 3/27

⎤ 0 ⎤ ⎡ 0 ⎥ ⎥ Tw 0 ⎥ 003 ⎥ ⎥⎢ ⎥ 0 ⎥ ⎢Tw ⎥ ⎢ 102 ⎥ 0 ⎥ ⎢Tw ⎥ ⎥⎢ ⎥ ⎥ 0 ⎥ ⎢ 201 ⎥ ⎢Tw ⎥ 0 ⎥ ⎢ 300 ⎥ ⎥ ⎥⎢ 0 ⎥ ⎢Tw ⎥ ⎥ ⎥ ⎢ 012 w 0 ⎥ ⎢T ⎥ ⎢ ⎥ 111 ⎥ 0 ⎥⎢ w ⎥ ⎢T ⎥ 0 ⎥ ⎥ ⎢ 210 ⎥ ⎢ w ⎥ 0 ⎥ ⎥ ⎢T021 ⎥ ⎥ ⎢ 1/8 ⎥ ⎥ ⎣Tw 120 ⎦ ⎥ 1/12⎥ Tw 030 1/18⎦ 1/27

(A.2)

98

J. Maestre et al. / Comput. Methods Appl. Mech. Engrg. 326 (2017) 70–101

w ˆw ezier patches Q and by permuting the indices Tiwjk to Tw ki j and T jki , the homogeneous coordinates for the two other B´ 2 w ˆ and Q3 can be obtained, respectively.

References [1] D.A. Andrews, P.S. Low, Role of red blood cells in thrombosis, Curr. Opin. Hematol. 6 (2) (1999) 76. [2] S.S. Bansode, S.K. Banarjee, D.D. Gaikwad, S.L. Jadhav, R.M. Thorat, Microencapsulation: a review, Int. J. Pharm. Sci. Rev. Res. 1 (2) (2010) 38–43. [3] D. Barthès-Biesel, Motion of a spherical microcapsule freely suspended in a linear shear flow, J. Fluid Mech. 100 (04) (1980) 831–853. [4] D. Barthès-Biesel, J.M. Rallison, The time-dependent deformation of a capsule freely suspended in a linear shear flow, J. Fluid Mech. 113 (1981) 251–267. [5] T. Gao, H.H. Hu, P.P. Castañeda, Rheology of a suspension of elastic particles in a viscous shear flow, J. Fluid Mech. 687 (2011) 209. [6] M.M. Villone, G. d’Avino, M.A. Hulsen, P.L. Maffettone, Dynamics of prolate spheroidal elastic particles in confined shear flow, Phys. Rev. E 92 (6) (2015) 062303. [7] K. Sugiyama, S. Ii, S. Takeuchi, S. Takagi, Y. Matsumoto, A full Eulerian finite difference approach for solving fluid–structure coupling problems, J. Comput. Phys. 230 (3) (2011) 596–627. [8] C.D. Eggleton, A.S. Popel, Large deformation of red blood cell ghosts in a simple shear flow, Phys. Fluids (1994-Present) 10 (8) (1998) 1834–1845. [9] P. Bagchi, P.C. Johnson, A.S. Popel, Computational fluid dynamic simulation of aggregation of deformable cells in a shear flow, J. Biomech. Eng. 127 (7) (2005) 1070–1080. [10] P. Bagchi, Mesoscale simulation of blood flow in small vessels, Biophys. J. 92 (6) (2007) 1858–1877. [11] C. Song, S.J. Shin, H.J. Sung, K.-S. Chang, Dynamic fluid-structure interaction of an elastic capsule in a viscous shear flow at moderate Reynolds number, J. Fluids Struct. 27 (3) (2011) 438–455. [12] W.-X. Huang, C.B. Chang, H.J. Sung, Three-dimensional simulation of elastic capsules in shear flow by the penalty immersed boundary method, J. Comput. Phys. 231 (8) (2012) 3340–3364. [13] Y. Sui, Y.-T. Chew, P. Roy, H.-T. Low, A hybrid method to study flow-induced deformation of three-dimensional capsules, J. Comput. Phys. 227 (12) (2008) 6351–6371. [14] R.M. MacMeccan, J.R. Clausen, G.P. Neitzel, C.K. Aidun, Simulating deformable particle suspensions using a coupled lattice-Boltzmann and finite-element method, J. Fluid Mech. 618 (2009) 13–39. [15] Y. Sui, X.B. Chen, Y.T. Chew, P. Roy, H.T. Low, Numerical simulation of capsule deformation in simple shear flow, Comput. & Fluids 39 (2) (2010) 242–250. [16] A. Kilimnik, W. Mao, A. Alexeev, Inertial migration of deformable capsules in channel flow, Phys. Fluids (1994-Present) 23 (12) (2011) 123302. [17] S.H. Schot, Eighty years of Sommerfeld’s radiation condition, His. Math. 19 (4) (1992) 385–401. [18] C. Pozrikidis, Finite deformation of liquid capsules enclosed by elastic membranes in simple shear flow, J. Fluid Mech. 297 (1995) 123–152. [19] S. Ramanujan, C. Pozrikidis, Deformation of liquid capsules enclosed by elastic membranes in simple shear flow: large deformations and the effect of fluid viscosities, J. Fluid Mech. 361 (1998) 117–143. [20] A. Leyrat-Maurin, D. Barthès-Biesel, Motion of a deformable capsule through a hyperbolic constriction, J. Fluid Mech. 279 (1994) 135–163. [21] C. Quéguiner, D. Barthès-Biesel, Axisymmetric motion of capsules through cylindrical channels, J. Fluid Mech. 348 (1997) 349–376. [22] A. Diaz, N. Pelekasis, D. Barthès-Biesel, Transient response of a capsule subjected to varying flow conditions: effect of internal fluid viscosity and membrane elasticity, Phys. Fluids (1994-Present) 12 (5) (2000) 948–957. [23] A. Diaz, D. Barthès-Biesel, Entrance of a bioartificial capsule in a pore, CMES Comput. Model. Eng. Sci. 3 (3) (2001) 321–338. [24] Y. Lefebvre, D. Barthès-Biesel, Motion of a capsule in a cylindrical tube: effect of membrane pre-stress, J. Fluid Mech. 589 (2007) 157–181. [25] E. Lac, D. Barthès-Biesel, N.A. Pelekasis, J. Tsamopoulos, Spherical capsules in three-dimensional unbounded Stokes flows: effect of the membrane constitutive law and onset of buckling, J. Fluid Mech. 516 (2004) 303–334. [26] E. Lac, D. Barthès-Biesel, Deformation of a capsule in simple shear flow: effect of membrane prestress, Phys. Fluids (1994-Present) 17 (7) (2005) 072105. [27] W.R. Dodson III, P. Dimitrakopoulos, Spindles, cusps, and bifurcation for capsules in Stokes flow, Phys. Rev. Lett. 101 (20) (2008) 208102. [28] W.R. Dodson, P. Dimitrakopoulos, Dynamics of strain-hardening and strain-softening capsules in strong planar extensional flows via an interfacial spectral boundary element algorithm for elastic membranes, J. Fluid Mech. 641 (2009) 263–296. [29] L. Zhu, J. Rabault, L. Brandt, The dynamics of a capsule in a wall-bounded oscillating shear flow, Phys. Fluids (1994-Present) 27 (7) (2015) 071902. [30] C. Rorai, A. Touchard, L. Zhu, L. Brandt, Motion of an elastic capsule in a constricted microchannel, Eur. Phys. J. E 38 (5) (2015) 1–13. [31] P.R. Zarda, S. Chien, R. Skalak, Interaction of viscous incompressible fluid with an elastic body, Comput. Methods Fluid-Solid Interact. Probl. (1977) 65–82. [32] R. Skalak, A. Tozeren, R.P. Zarda, S. Chien, Strain energy function of red blood cell membranes, Biophys. J. 13 (3) (1973) 245. [33] R. Skalak, N. Ozkaya, T.C. Skalak, Biofluid mechanics, Annu. Rev. Fluid Mech. 21 (1) (1989) 167–200. [34] G. Pieper, H. Rehage, D. Barthès-Biesel, Deformation of a capsule in a spinning drop apparatus, J. Colloid Interface Sci. 202 (2) (1998) 293–300. [35] M. Carin, D. Barthès-Biesel, F. Edwards-Lévy, C. Postel, D.C. Andrei, Compression of biocompatible liquid-filled HSA-alginate capsules: Determination of the membrane mechanical properties, Biotechnol. Bioeng. 82 (2) (2003) 207–212.

J. Maestre et al. / Comput. Methods Appl. Mech. Engrg. 326 (2017) 70–101

99

[36] M. Husmann, H. Rehage, E. Dhenin, D. Barthès-Biesel, Deformation and bursting of nonspherical polysiloxane microcapsules in a spinningdrop apparatus, J. Colloid Interface Sci. 282 (1) (2005) 109–119. [37] C. Pozrikidis, Numerical simulation of the flow-induced deformation of red blood cells, Ann. Biomed. Eng. 31 (10) (2003) 1194–1205. [38] H. Zhao, A.H. Isfahani, L.N. Olson, J.B. Freund, A spectral boundary integral method for flowing blood cells, J. Comput. Phys. 229 (10) (2010) 3726–3744. [39] Z. Peng, R.J. Asaro, Q. Zhu, Multiscale simulation of erythrocyte membranes, Phys. Rev. E 81 (2010) 031904. [40] Z. Peng, Y.-L. Chen, H. Lu, Z. Pan, H.-C. Chang, Mesoscale simulations of two model systems in biophysics: from red blood cells to DNAs, Comput. Part. Mech. 2 (4) (2015) 339–357. [41] C. Dupont, A.-V. Salsac, D. Barthès-Biesel, M. Vidrascu, P. Le Tallec, Influence of bending resistance on the dynamics of a spherical capsule in shear flow, Phys. Fluids (1994-Present) 27 (5) (2015) 051902. [42] D. Barthès-Biesel, Motion and deformation of elastic capsules and vesicles in flow, Annu. Rev. Fluid Mech. 48 (2016) 25–52. [43] J. Walter, A.-V. Salsac, D. Barthès-Biesel, P. Le Tallec, Coupling of finite element and boundary integral methods for a capsule in a Stokes flow, Internat. J. Numer. Methods Engrg. 83 (7) (2010) 829–850. [44] E. Foessel, J. Walter, A.-V. Salsac, D. Barthès-Biesel, Influence of internal viscosity on the large deformation and buckling of a spherical capsule in a simple shear flow, J. Fluid Mech. 672 (2011) 477–486. [45] J. Walter, A.-V. Salsac, D. Barthès-Biesel, Ellipsoidal capsules in simple shear flow: prolate versus oblate initial shapes, J. Fluid Mech. 676 (2011) 318–347. [46] X.-Q. Hu, A.-V. Salsac, D. Barthès-Biesel, Flow of a spherical capsule in a pore with circular or square cross-section, J. Fluid Mech. 705 (2012) 176–194. [47] T. Omori, T. Ishikawa, Y. Imai, T. Yamaguchi, Membrane tension of red blood cells pairwisely interacting in simple shear flow, J. Biomech. 46 (3) (2013) 548–553. [48] C. Dupont, F. Delahaye, A.-V. Salsac, D. Barthès-Biesel, Off-plane motion of an oblate capsule in a simple shear flow, Comput. Methods Biomech. Biomed. Eng. 16 (Suppl. 1) (2013) 4–5. [49] T.J.R. Hughes, J.A. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Comput. Methods Appl. Mech. Engrg. 194 (39–41) (2005) 4135–4195. [50] Y. Bazilevs, V.M. Calo, Y. Zhang, T.J.R. Hughes, Isogeometric fluid–structure interaction analysis with applications to arterial blood flow, Comput. Mech. 38 (4–5) (2006) 310–322. [51] Y. Zhang, Y. Bazilevs, S. Goswami, C.L. Bajaj, T.J.R. Hughes, Patient-specific vascular NURBS modeling for isogeometric analysis of blood flow, Comput. Methods Appl. Mech. Engrg. 196 (29) (2007) 2943–2959. [52] Y. Bazilevs, J.R. Gohean, T.J.R. Hughes, R.D. Moser, Y. Zhang, Patient-specific isogeometric fluid–structure interaction analysis of thoracic aortic blood flow due to implantation of the Jarvik 2000 left ventricular assist device, Comput. Methods Appl. Mech. Engrg. 198 (45) (2009) 3534–3550. [53] V. Chivukula, J. Mousel, J. Lu, S. Vigmostad, Micro-scale blood particulate dynamics using a non-uniform rational B-spline-based isogeometric analysis, Int. J. Numer. Methods Biomed. Eng. 30 (12) (2014) 1437–1459. [54] C. Politis, A. Ginnis, P. Kaklis, K. Belibassakis, C. Feurer, An isogeometric BEM for exterior potential-flow problems in the plane, Joint Conf. Geom. Phys. Model. (2009) 349–354. [55] R.N. Simpson, S.P.A. Bordas, J. Trevelyan, T. Rabczuk, A two-dimensional isogeometric boundary element method for elastostatic analysis, Comput. Methods Appl. Mech. Engrg. 209–212 (2012) 87–100. [56] K. Li, X. Qian, Isogeometric analysis and shape optimization via boundary integral, Comput. Aided Des. 43 (11) (2011) 1427–1437. [57] K. Belibassakis, T. Gerostathis, K. Kostas, C. Politis, P. Kaklis, A. Ginnis, C. Feurer, A BEM-isogeometric method for the ship waveresistance problem, Ocean Eng. 60 (2013) 53–67. [58] M.J. Peake, J. Trevelyan, G. Coates, Extended isogeometric boundary element method (XIBEM) for two-dimensional Helmholtz problems, Comput. Methods Appl. Mech. Engrg. 259 (2013) 93–102. [59] C. Politis, A. Papagiannopoulos, K. Belibassakis, P. Kaklis, K. Kostas, A. Ginnis, T. Gerostathis, An isogeometric BEM for exterior potentialflow problems around lifting bodies, in: 11th World Congress on Computational Mechanics, WCCM XI, International Center for Numerical Methods in Engineering (CIMNE), 2014, pp. 2433–2444. [60] M.J. Peake, J. Trevelyan, G. Coates, Extended isogeometric boundary element method (XIBEM) for three-dimensional medium-wave acoustic scattering problems, Comput. Methods Appl. Mech. Engrg. 284 (2015) 762–780. [61] A. Aimi, M. Diligenti, M. Sampoli, A. Sestini, Isogemetric analysis and symmetric Galerkin BEM: A 2D numerical study, Appl. Math. Comput. 272 (2016) 173–186. [62] A. Joneidi, C. Verhoosel, P. Anderson, Isogeometric boundary integral analysis of drops and inextensible membranes in isoviscous flow, Comput. & Fluids 109 (2015) 49–66. [63] L. Heltai, M. Arroyo, A. DeSimone, Nonsingular isogeometric boundary element method for Stokes flows in 3D, Comput. Methods Appl. Mech. Engrg. 268 (2014) 514–539. [64] L. Heltai, J. Kiendl, A. DeSimone, A. Reali, A natural framework for isogeometric fluid–structure interaction based on BEM–shell coupling, Comput. Methods Appl. Mech. Engrg. 316 (2017) 522–546. [65] T.W. Sederberg, J. Zheng, A. Bakenov, A. Nasri, T-splines and T-NURCCs, ACM Trans. Graph. 22 (3) (2003) 477–484. [66] M.A. Scott, X. Li, T.W. Sederberg, T.J.R. Hughes, Local refinement of analysis-suitable T-splines, Comput. Methods Appl. Mech. Engrg. 213 (2012) 206–222. [67] X. Li, J. Zheng, T.W. Sederberg, T.J.R. Hughes, M.A. Scott, On linear independence of T-spline blending functions, Comput. Aided Geom. Design 29 (1) (2012) 63–76.

100

J. Maestre et al. / Comput. Methods Appl. Mech. Engrg. 326 (2017) 70–101

[68] M.J. Borden, M.A. Scott, J.A. Evans, T.J.R. Hughes, Isogeometric finite element data structures based on Bézier extraction of NURBS, Internat. J. Numer. Methods Engrg. 87 (1–5) (2011) 15–47. [69] M.A. Scott, M.J. Borden, C.V. Verhoosel, T.W. Sederberg, T.J.R. Hughes, Isogeometric finite element data structures based on Bézier extraction of T-splines, Internat. J. Numer. Methods Engrg. 88 (2) (2011) 126–156. [70] M.A. Scott, R.N. Simpson, J.A. Evans, S. Lipton, S.P.A. Bordas, T.J.R. Hughes, T.W. Sederberg, Isogeometric boundary element analysis using unstructured T-splines, Comput. Methods Appl. Mech. Engrg. 254 (0) (2013) 197–221. [71] R. Dimitri, L. De Lorenzis, M.A. Scott, P. Wriggers, R.L. Taylor, G. Zavarise, Isogeometric large deformation frictionless contact using T-splines, Comput. Methods Appl. Mech. Engrg. 269 (2014) 394–414. [72] R.N. Simpson, M.A. Scott, M. Taus, D.C. Thomas, H. Lian, Acoustic isogeometric boundary element analysis, Comput. Methods Appl. Mech. Engrg. 269 (2014) 265–290. [73] K.V. Kostas, A.I. Ginnis, C.G. Politis, P.D. Kaklis, Ship-hull shape optimization with a T-spline based BEM-isogeometric solver, Comput. Methods Appl. Mech. Engrg. 284 (2015) 611–622. [74] J. Maestre, I. Cuesta, J. Pallares, An unsteady 3d isogeometrical boundary element analysis applied to nonlinear gravity waves, Comput. Methods Appl. Mech. Engrg. 310 (2016) 112–133. [75] Y. Bazilevs, V.M. Calo, J.A. Cottrell, J.A. Evans, T.J.R. Hughes, S. Lipton, M.A. Scott, T.W. Sederberg, Isogeometric analysis using Tsplines, Comput. Methods Appl. Mech. Engrg. 199 (5–8) (2010) 229–263. [76] D. Schillinger, L. Dedè, M.A. Scott, J.A. Evans, M.J. Borden, E. Rank, T.J.R. Hughes, An isogeometric design-through-analysis methodology based on adaptive hierarchical refinement of NURBS, immersed boundary methods, and T-spline CAD surfaces, Comput. Methods Appl. Mech. Engrg. 249–252 (2012) 116–150. [77] D.J. Benson, Y. Bazilevs, E. de Luycker, M.C. Hsu, M.A. Scott, T.J.R. Hughes, T. Belytschko, A generalized finite element formulation for arbitrary basis functions: From isogeometric analysis to XFEM, Internat. J. Numer. Methods Engrg. 83 (6) (2010) 765–785. [78] C.V. Verhoosel, M.A. Scott, R. De Borst, T.J.R. Hughes, An isogeometric approach to cohesive zone modeling, Internat. J. Numer. Methods Engrg. 87 (1–5) (2011) 336–360. [79] C.V. Verhoosel, M.A. Scott, T.J.R. Hughes, R. de Borst, An isogeometric analysis approach to gradient damage models, Internat. J. Numer. Methods Engrg. 86 (1) (2011) 115–134. [80] M.J. Borden, C.V. Verhoosel, M.A. Scott, T.J.R. Hughes, C.M. Landis, A phase-field description of dynamic brittle fracture, Comput. Methods Appl. Mech. Engrg. 217–220 (2012) 77–95. [81] S.S. Ghorashi, N. Valizadeh, S. Mohammadi, T. Rabczuk, T-spline based XIGA for fracture analysis of orthotropic media, Comput. Struct. 147 (2015) 138–146. [82] Y. Bazilevs, M.C. Hsu, M.A. Scott, Isogeometric fluid-structure interaction analysis with emphasis on non-matching discretizations, and with application to wind turbines, Comput. Methods Appl. Mech. Engrg. 249–252 (2012) 28–41. [83] M.-C. Hsu, D. Kamensky, F. Xu, J. Kiendl, C. Wang, M. Wu, J. Mineroff, A. Reali, Y. Bazilevs, M. Sacks, Dynamic and fluidstructure interaction simulations of bioprosthetic heart valves using parametric design with T-splines and Fung-type material models, Comput. Mech. 55 (6) (2015) 1211–1225. [84] A. Buffa, G. Sangalli, R. Vázquez, Isogeometric methods for computational electromagnetics: B-spline and T-spline discretizations, J. Comput. Phys. 257 (PB) (2014) 1291–1320. [85] G. Lorenzo, M.A. Scott, K. Tew, T.J.R. Hughes, H. Gomez, Hierarchically refined and coarsened splines for moving interface problems, with particular application to phase-field models of prostate tumor growth, Comput. Methods Appl. Mech. Engrg. (2017). [86] J. Deng, F. Chen, X. Li, C. Hu, W. Tong, Z. Yang, Y. Feng, Polynomial splines over hierarchical T-meshes, Graph. Models 70 (4) (2008) 76–86. [87] M.A. Scott, D.C. Thomas, E.J. Evans, Isogeometric spline forests, Comput. Methods Appl. Mech. Engrg. 269 (2014) 222–264. [88] E.J. Evans, M.A. Scott, X. Li, D.C. Thomas, Hierarchical T-splines: Analysis-suitability, Bézier extraction, and application as an adaptive basis for isogeometric analysis, Comput. Methods Appl. Mech. Engrg. 284 (2015) 1–20. [89] C. Giannelli, B. Jüttler, H. Speleers, THB-splines: The truncated basis for hierarchical splines, Comput. Aided Geom. Design 29 (7) (2012) 485–498. [90] P. Hennig, S. Müller, M. Kästner, Bézier extraction and adaptive refinement of truncated hierarchical NURBS, Comput. Methods Appl. Mech. Engrg. 305 (2016) 316–339. [91] X. Wei, Y. Zhang, T.J.R. Hughes, M.A. Scott, Truncated hierarchical Catmull–Clark subdivision with local refinement, Comput. Methods Appl. Mech. Engrg. 291 (2015) 1–20. [92] S.K. Boey, D.H. Boal, D.E. Discher, Simulations of the erythrocyte cytoskeleton at large deformation. I. Microscopic models, Biophys. J. 75 (3) (1998) 1573–1583. [93] D.E. Discher, D.H. Boal, S.K. Boey, Simulations of the erythrocyte cytoskeleton at large deformation. II. Micropipette aspiration, Biophys. J. 75 (3) (1998) 1584–1597. [94] C. Pozrikidis, Numerical simulation of the flow-induced deformation of red blood cells, Ann. Biomed. Eng. 31 (10) (2003) 1194–1205. [95] C. Pozrikidis, Numerical simulation of cell motion in tube flow, Ann. Biomed. Eng. 33 (2) (2005) 165–178. [96] I.V. Pivkin, G.E. Karniadakis, Accurate coarse-grained modeling of red blood cells, Phys. Rev. Lett. 101 (11) (2008) 118105. [97] D.A. Fedosov, B. Caswell, G.E. Karniadakis, A multiscale red blood cell model with accurate mechanics, rheology, and dynamics, Biophys. J. 98 (10) (2010) 2215–2225. [98] H. Zhao, A.H. Isfahani, L.N. Olson, J.B. Freund, A spectral boundary integral method for flowing blood cells, J. Comput. Phys. 229 (10) (2010) 3726–3744. [99] X. Li, Z. Peng, H. Lei, M. Dao, G.E. Karniadakis, Probing red blood cell mechanics, rheology and dynamics with a two-component multiscale model, Phil. Trans. R. Soc. A 372 (2021) (2014) 20130389.

J. Maestre et al. / Comput. Methods Appl. Mech. Engrg. 326 (2017) 70–101

101

[100] D. Barthès-Biesel, A. Diaz, E. Dhenin, Effect of constitutive laws for two-dimensional membranes on flow-induced capsule deformation, J. Fluid Mech. 460 (2002) 211–222. [101] X. Li, M.A. Scott, Analysis-suitable T-splines: characterization, refineability, and approximation, Math. Models Methods Appl. Sci. 24 (06) (2014) 1141–1164. [102] M. Taus, G.J. Rodin, T.J. Hughes, Isogeometric analysis of boundary integral equations: High-order collocation methods for the singular and hyper-singular equationsac, Math. Models Methods Appl. Sci. 26 (08) (2016) 1447–1480. [103] T. Takacs, B. Jüttler, Existence of stiffness matrix integrals for singularly parameterized domains in isogeometric analysis, Comput. Methods Appl. Mech. Engrg. 200 (49) (2011) 3568–3582. [104] S.-M. Hu, Conversion between triangular and rectangular Bézier patches, Comput. Aided Geom. Design 18 (7) (2001) 667–671. [105] L. Yan, X. Han, J. Liang, Conversion between triangular Bézier patches and rectangular Bézier patches, Appl. Math. Comput. 232 (2014) 469–478. [106] C. Pozrikidis, A Practical Guide to Boundary Element Methods with the Software Library BEMLIB, CRC Press, 2002. [107] V. Mantic, A new formula for the C-matrix in the Somigliana identity, J. Elasticity 33 (3) (1993) 191–201. [108] C.A. Brebbia, J. Domínguez, Boundary Elements: An Introductory Course, second ed., WIT Press, 1992. Computational Mechanics. [109] J.C. Lachat, J.O. Watson, Effective numerical treatment of boundary integral equations: A formulation for three-dimensional elastostatics, Internat. J. Numer. Methods Engrg. 10 (5) (1976) 991–1005. [110] J. Dominguez, Boundary Elements in Dynamics, Wit Press, 1993. [111] F. Auricchio, F. Calabrò, T. Hughes, A. Reali, G. Sangalli, A simple algorithm for obtaining nearly optimal quadrature rules for NURBSbased isogeometric analysis, Comput. Methods Appl. Mech. Engrg. 249–252 (2012) 15–27. [112] T. Hughes, A. Reali, G. Sangalli, Efficient quadrature for NURBS-based isogeometric analysis, Comput. Methods Appl. Mech. Engrg. 199 (5–8) (2010) 301–313. [113] T. Cruse, R. Aithal, Non-singular boundary integral equation implementation, Internat. J. Numer. Methods Engrg. 36 (2) (1993) 237–254. [114] Q. Huang, T. Cruse, Some notes on singular integral techniques in boundary element analysis, Internat. J. Numer. Methods Engrg. 36 (15) (1993) 2643–2659. [115] P.R. Johnston, D. Elliott, A sinh transformation for evaluating nearly singular boundary element integrals, Internat. J. Numer. Methods Engrg. 62 (4) (2005) 564–578. [116] J. Lekner, Viscous flow through pipes of various cross-sections, Eur. J. Phys. 28 (3) (2007) 521. [117] T. Sederberg, D. Cardon, J. Zheng, T. Lyche, T-spline simplification and local refinement, in: SIGGRAPH 2004, ACM, 2004, pp. 276–283. [118] D.C. Thomas, M.A. Scott, J.A. Evans, K. Tew, E.J. Evans, Bézier projection: A unified approach for local projection and quadrature-free refinement and coarsening of NURBS and T-splines with particular application to isogeometric design and analysis, Comput. Methods Appl. Mech. Engrg. 284 (2015) 55–105. Isogeometric Analysis Special Issue. [119] C. Pozrikidis, Effect of membrane bending stiffness on the deformation of capsules in simple shear flow, J. Fluid Mech. 440 (2001) 269–291. [120] Y. Wang, P. Dimitrakopoulos, A three-dimensional spectral boundary element algorithm for interfacial dynamics in Stokes flow, Phys. Fluids (1994-Present) 18 (8) (2006) 082106. [121] A.T. Chwang, T.Y.-T. Wu, Hydromechanics of low-Reynolds-number flow. Part 2. Singularity method for Stokes flows, J. Fluid Mech. 67 (04) (1975) 787–815. [122] B.J. Bentley, L.G. Leal, An experimental investigation of drop deformation and breakup in steady, two-dimensional linear flows, J. Fluid Mech. 167 (1986) 241–283. [123] S.M. Avtan, M. Kaya, N. Orhan, A. Arslan, N. Arican, A.S. Toklu, C. Gürses, I. Elmas, M. Kucuk, B. Ahishali, The effects of hyperbaric oxygen therapy on blood–brain barrier permeability in septic rats, Brain Res. 1412 (2011) 63–72. [124] H.L. Goldsmith, Red cell motions and wall interactions in tube flow, in: Federation Proceedings, vol. 30, 1971, pp. 1578–1590. [125] A. Helmy, D. Barthès-Biesel, Migration of a spherical capsule freely suspended in an unbounded parabolic flow, J. Mec. Theor. Appl. 1 (5) (1982) 859–880. [126] G. Coupier, B. Kaoui, T. Podgorski, C. Misbah, Noninertial lateral migration of vesicles in bounded Poiseuille flow, Phys. Fluids (1994Present) 20 (11) (2008) 111702. [127] S.K. Doddi, P. Bagchi, Lateral migration of a capsule in a plane Poiseuille flow in a channel, Int. J. Multiph. Flow. 34 (10) (2008) 966–986. [128] L. Shi, T.-W. Pan, R. Glowinski, Numerical simulation of lateral migration of red blood cells in Poiseuille flows, Internat. J. Numer. Methods Fluids 68 (11) (2012) 1393–1408. [129] R.K. Singh, X. Li, K. Sarkar, Lateral migration of a capsule in plane shear near a wall, J. Fluid Mech. 739 (2014) 421–443. [130] S.-M. Hu, Conversion of a triangular Bézier patch into three rectangular Bézier patches, Comput. Aided Geom. Design 13 (3) (1996) 219–226.