A partition of unity approach to fluid mechanics and fluid-structure interaction
Maximilian Balmus, Andre Massing, Johan Hoffman, Reza Razavi, David Nordsletten
AA partition of unity approach to fluid mechanics andfluid-structure interaction
Maximilian Balmus a, ∗ , Andr´e Massing b,c , Johan Hoffman d , Reza Razavi a ,David A. Nordsletten a,e a Department of Biomedical Engineering, School of Imaging Sciences and Biomedical Engineering,King’s College London, King’s Health Partners, London SE1 7EH, United Kingdom b Department of Mathematics and Mathematical Statistics, Ume˚a University, SE-90187 Ume˚a, Sweden c Department of Mathematical Sciences, Norwegian University of Science and Technology, NO-7491Trondheim, Norway d Division of Computational Science and Technology, KTH Royal Institute of Technology, Sweden e Department of Biomedical Engineering and Cardiac Surgery, University of Michigan, MI, USA
Abstract
For problems involving large deformations of thin structures, simulating fluid-structureinteraction (FSI) remains challenging largely due to the need to balance computationalfeasibility, efficiency, and solution accuracy. Overlapping domain techniques have beenintroduced as a way to combine the fluid-solid mesh conformity, seen in moving-meshmethods, without the need for mesh smoothing or re-meshing, which is a core character-istic of fixed mesh approaches. In this work, we introduce a novel overlapping domainmethod based on a partition of unity approach. Unified function spaces are defined as aweighted sum of fields given on two overlapping meshes. The method is shown to achieveoptimal convergence rates and to be stable for steady-state Stokes, Navier-Stokes, andALE Navier-Stokes problems. Finally, we present results for FSI in the case of a 2Dmock aortic valve simulation. These initial results point to the potential applicability ofthe method to a wide range of FSI applications, enabling boundary layer refinement andlarge deformations without the need for re-meshing or user-defined stabilization.
Keywords:
Finite element methods, Fluid-structure interaction, Overlapping domains,Partition of unity ∗ Corresponding author
Email addresses: [email protected] (Maximilian Balmus ), [email protected] (Andr´e Massing), [email protected] (Johan Hoffman), [email protected] (David A. Nordsletten)
Preprint submitted to Computer Methods in Applied Mechanics and Engineering February 19, 2019 a r X i v : . [ c s . C E ] F e b . Introduction Fluid-structure interaction (FSI) problems involving thin solids which undergo largedeformations and translations can be encountered in a significant number of engineeringapplications. In industry, we have examples such as the design of parachutes [1, 2, 3]and wind-turbines [4, 5]. In cardiovascular research, the simulation of valves [6, 7] andimplanted devices [8] holds a great potential for better understanding and treatment ofa number of pathologies such as valve stenosis, regurgitation, heart failure and outflowobstruction. Building such models, however, remains challenging due to the need tobalance computational costs and solution accuracy. In the case of cardiac valves, forexample, studies typically require simplifications of the domain’s geometry [9] or thefluid models [10].FSI approaches for this class of problems can be grouped into three main categories:interface-tracking, interface-capturing and overlapping domains methods [11, 12]. Inthe case of interface-tracking, the fluid problem is typically based on the ArbitraryLagrangian-Eulerian (ALE) [13, 14, 15] formulation. This allows for the fluid domainto deform with the solid and enables adjusting the element resolution close to surfacesin order to more accurately represent boundary layers. In [16], an ALE based methodis shown to produce superior results when considering moderate deformations to threevariations of fictitious domain (examples of interface-capturing) for similar mesh reso-lutions. However, for large deformations, it is known that the distortion of the fluidmesh can diminish the quality of elements and negatively impact the accuracy of solu-tions [11]. While multiple re-meshing techniques [17] have been proposed, the processintroduces grid interpolation errors and its computational cost effectiveness is linked tothe frequency at which the mesh needs to be adjusted. As shown in [18], where theparticular case of cardiac valves is discussed, the ALE based simulations take more timeto run than competing interface-capturing methods.In contrast, interface-capturing methods do not require boundary fitted meshes forthe fluid and solid and thus avoid the need for re-meshing. Examples include the Fic-titious Domain Method [19, 20] (FDM), where the coupling between the solid and fluidis achieved via additional Lagrange multiplier terms [21], and the Immersed Finite El-ement method (IFEM) [22, 23, 24, 25], where the kinematic constraints between thefluid and solid is imposed through interpolation and distribution of local body forces.While both FDM and IFEM require the construction of a solid mesh, in the ImmersedStructural Potential Method (ISPM) [26] both problems are solved on the same mesh,with the solid being represented as a moving collection of quadrature points. However,these methods lack a conforming interface between the fluid and solid and can resultin poor approximations of pressure jumps and surface stresses [27, 28]. For this reason,mesh-adaptation [29] and XFEM enrichment [30, 31, 32, 33, 34, 35, 36, 37] techniqueshave been proposed in order to overcome these issues. Also a significant challenge forinterface-capturing, as noted in [38], remains the fact the resolution of the fluid flowaround the structural surface is limited by the local element size of the fluid mesh. Inpractice, this leads to over-refinement of the fluid mesh along the moving trajectory ofthe solid and to significant increases of the computational costs.Overlapping domain techniques [39, 40, 41, 42, 43, 44, 45, 46, 12, 47] have beenproposed with the aim of combining the advantages of interface-tracking and interfacecapturing techniques: fluid-solid mesh conformity, boundary layer tracking and eliminat-2ng re-meshing. The crux of these methods is the decomposition of the fluid probleminto a background coarse component and solution-enriching embedded component thatenvelops the structures. A challenging aspect, however, is the coupling of the two fluiddomains which has to be done weakly due the non-matching fluid-fluid interface. Theuse of Lagrange multipliers for example, see [30], is impeded by the need to properlychoose function spaces in order to guarantee that the inf-sup condition holds for ar-bitrary moving interfaces. Alternatively, stabilization techniques can be employed tocircumvent the inf-sup condition [48, 49, 50]. Overlapping mesh methods for the Stokesproblem [51, 52, 46] which use Nitsche’s method [53] avoid introducing an additionalLagrange multiplier field, but nevertheless they require additional, parameter-dependentstabilization terms for the velocity and pressure jump in the vicinity of the fluid-fluid in-terface to guarantee optimal convergence rates and good system conditioning irrespectiveof the particular overlap configuration.In this paper, we propose a new flexible and robust overlapping domain method whichuses the partition of unity (PUFEM) approach [54, 55] to decompose the fluid domain intoa background mesh and an embedded mesh which can overlap in an arbitrary manner.On each mesh, standard mixed and inf-sup stable velocity-pressure function spaces usingTaylor-Hood elements are defined. In the final finite element formulation of the fluidproblem, a unified global function space is then used which is defined by taking a properlyweighted sum of the function spaces associated with each mesh. To avoid ill-conditioningof the resulting system, additional constraints are introduced in (parts of) the overlapregion.The rest of this paper is structured as follows. After briefly recalling the classi-cal mixed finite element approach for the Stokes problem in Section 2.1, we introduceits PUFEM based overlapping domain formulation including a detailed description ofthe domain set-up, weighted function spaces, and imposed constraints, see Section 2.2,followed by a short discussion of some computational aspects in Section 2.3. Then inSection 3, we explain how to combine the PUFEM approach with an ALE formulation ofthe Navier-Stokes problem to treat moving fluid domain and fluid–structure interactionproblems. In Section 4, a number of numerical experiments are conducted. First, weinvestigate the stability and accuracy of our PU approach for a number of fluid flowproblems posed on fixed static domains, see Section 4.1–4.3. To demonstrate the capa-bility to handle large changes in the fluid domain geometry, we then consider a fluid flowdriven by a oscillating cylinder in Section 4.4 before we turn to a full FSI problem inSection 4.5, where we compare a classical ALE-based approach with our novel combinedALE-PUFEM discretization for a two-dimensional mock aortic valve simulation. Finallyin Section 5, we summarize our results and discuss potential future developments.
2. A Partition of unity finite element method for the Stokes problem
In this section, we introduce the main concepts and the basic setup for PUFEM. Webegin by reviewing the classic mixed FEM approach to Stokes flow (Section 2.1) andpresent the key differences introduced in the PUFEM approach (Section 2.2). Finally,in Section 2.3, we describe the process by which we identify the polygonal intersectionsof overlapping elements and perform the necessary integrations.3 .1. Classical mixed FEM approach to Stokes problems
Let Ω f ⊂ R d be an arbitrary fluid domain on which we solve our problem. Since ourfocus is oriented towards FSI simulations, we also introduce a solid domain Ω s whichfor now, we consider to be fixed and rigid. The boundary of the fluid problem, Γ, iscomposed of three, non-overlapping regions: the portions where Dirichlet and Neumannboundary conditions are applied, (Γ D and Γ N , respectively), and the fluid-solid interface,Γ fs . We then can write Stokes problem as: find the ( v , p ) such that, µ ∇ v − ∇ p = in Ω f , (1a) ∇ · v = 0 in Ω f , (1b)( µ ∇ v − p I ) · n = t N on Γ N , (1c) v = v D on Γ D , (1d) v = on Γ fs , (1e)where µ is the fluid dynamic viscosity constant. For simplicity, we consider that theNeumann boundary condition tractions are null, i.e. t N = .In the case of the Stokes problem in (1), the classic continuous weak form can bewritten as: find ( v , p ) ∈ V D × W such that, (cid:90) Ω µ ∇ v : ∇ w − p ∇ · w + q ∇ · v d Ω = 0 , ∀ ( w , q ) ∈ V × W , (2)where V D , V ⊂ H (Ω f ) and W = L (Ω f ). Here the subscripts D and 0 indicatethat the V D , V ⊂ V subspaces are built such that they incorporate the Dirichlet andzero boundary value conditions on Γ D . Proofs of the well-posedness of this problemcan be found in [56] and [57]. In the discrete setting, the resulting weak form is: find( v h , p h ) ∈ V hD × W h such that, (cid:90) Ω hf µ ∇ v h : ∇ w h − p h ∇ · w h + q h ∇ · v h d Ω = 0 , ∀ ( w h , q h ) ∈ V h × W h . (3)In this study, we will use the classic approach in (3) to compare with the PUFEMapproach. In this case, we use the LBB stable P − P Taylor-Hood elements [58]. Thus,the discrete test and trial function spaces can be defined as: V h = { v h ∈ C (Ω hf ) | v | τ ∈ (cid:2) P ( τ ) (cid:3) d for τ ∈ T f } , (4a) W h = { p h ∈ C (Ω hf ) | p | τ ∈ P ( τ ) for τ ∈ T f } . (4b)Here Ω hf is a discrete equivalent of Ω f , defined by a tessellation T f . h denotes theelement size defined as the diameter of the circumcircle. In preparation for the PUFEMdiscussion, we can also define Ω hs as the discrete solid domain. P k designates the set ofpolynomial function of order k .Previous works [56, 57] derive a priori error estimates where: || v − v h || + || p − p h || ≤ C inf w h ∈ V h || v − w h || + C inf q h ∈W h || p − q h || , (5)4here C and C are positive constants independent of h , and || · || and || · || denotethe L and H norms, respectively. Using P − P elements, if ( v , p ) ∈ H (Ω) × H (Ω),then from interpolation theory [56]: || v − v h || + || p − p h || ≤ Ch ( | v | + | p | ) . (6) The core difference between the classic approach and the PUFEM setup is that thelatter is composed of two overlapping meshes: background and embedded (see Fig. 1),and both have a corresponding discrete domain over which they are defined. Thus, weassume the background domain, Ω hb = Ω hf ∪ Ω hs , encompasses the entirety of the discretefluid and solid domains. The embedded domain, Ω he , which satisfies Ω hs ⊂ Ω he ⊆ Ω hb , isdesigned to incorporate the solid and extends into fluid domains providing a boundarylayer. The fluid and solid regions of Ω he are separated by the Γ h fs interface. Additionally,Γ h ff is the outer boundary of Ω he and serves as its interface with Ω hb . The two mesheseach have a corresponding tessellation, T b and T e , and element sizes denoted by h b and h e , respectively. Let N vb and N pb be the set of nodes defining the background mesh in thecase of quadratic and linear interpolations. N ve and N pe are their embedded analogues.For the remainder of this paper, we consider the case where the meshes are 2D andtriangular.The discrete weak form for the PUFEM approach can be written as: find ( v h , p h ) ∈ V hD × W h such that, (cid:90) Ω hb µ ∇ v h : ∇ w h − p h ∇ · w h + q h ∇ · v h d Ω = 0 , ∀ ( w h , q h ) ∈ V h × W h (7)where V hD , V h and W h are PUFEM function spaces that represent the counterparts of V hD , V h and W h introduced in Section 2.1. The PUFEM function spaces are definedas the weighted sums of function spaces with support on the background and embeddedmeshes. Thus, based on the same notation of Melenk and Babuska [54], the spaces canbe written as: V h = (cid:0) − ψ h (cid:1) V hb, ∗ + ψ h V he, ∗ (8a) W h = (cid:0) − ψ h (cid:1) W hb, ∗ + ψ h W he, ∗ (8b)where the local spaces are defined as: V hk, ∗ = { v hk ∈ V hk | v k | x = v m | x for x ∈ X vk } (9a) W hk, ∗ = { p hk ∈ W hk | p k | x = p m | x for x ∈ X pk } (9b)for k, m ∈ { b, e } , k (cid:54) = m , and V hk = { v h ∈ C (Ω hk ) | v h | τ ∈ (cid:2) P ( τ ) (cid:3) d ∀ τ ∈ T k } for k ∈ { b, e } (10a) W hb = { p h ∈ C (Ω hg ) | p h | τ ∈ P ( τ ) ∀ τ ∈ T b } (10b) W he = { p h ∈ L (Ω he ) ∩ (cid:2) C (Ω hs ) ⊕ C (Ω he ∩ Ω hf ) (cid:3) | p h | τ ∈ P ( τ ) ∀ τ ∈ T e } (10c)5ote that with the exception of W he , all the discrete local spaces are continuous.We allow the functions in W he to be discontinuous across Γ h fs in order to be able torepresent the pressure jump between the fluid and embedded solid. The X vk and X pk sets contain nodes that are constrained in order to avoid ill-conditioning. More detailson how these sets are constructed can be found in Section 2.2.1.We define the field ψ h ∈ V he ( V he is the scalar equivalent of V he ) such that supp( ψ h ) ⊂ Ω he , its codomain is [0 ,
1] and ψ h = 0 on Γ h ff . For notational simplicity, we considerthat all fields defined on the embedded mesh (i.e. ψ h , v he and p he ) are zero outside theregion of overlap. In practice, we built ψ h using the following series of steps. First,we solved a diffusion problem on the embedded mesh, where we constrain the field tobe zero on Γ h ff and ten on Γ h fs . Subsequently, we apply a upper boundary thresholdsuch that the maximum node value is one. Finally, we smoothen the field using thehermitian polynomial f ( u ) = − u + 3 u . An example of the resulting field in the caseof circular mesh can be seen Fig. 2. This approach was used irrespective of element sizefor consistency. Note however, that this process was chosen on an ad hoc basis. A fullinvestigation on the effects of the support area and shape of ψ remain to be investigated.One of the benefits that we derive from the weighted form of the PUFEM is that itguarantees a smooth transition across Γ h ff , irrespective of the mesh size used in eitherbackground or embedded mesh. As a result, this should always prevent having mis-matched fluxes (i.e. when computed on both sides of the boundary) and the loss of massacross the interface. Previously we introduced four sets of nodes that we want to constrain X vb , X ve , X pb and X pe and in this subsection we will discuss both their necessity and how we decidewhich nodes belong to them. With the introduction of the weighting field, the standardbasis function support concept becomes insufficient to understand the contribution ofeach DOF to the total solution. Thus, we introduce a new metric of effective supportfraction which for an embedded DOF we define as: E ( x i ) = (cid:34)(cid:90) Ω hf ( ψ h ˆ φ i ) d Ω (cid:35) (cid:34)(cid:90) Ω hf ( ˆ φ i ) d Ω (cid:35) − (11)where ˆ φ i is the basis function of a degree of freedom corresponding to the x i node co-ordinate. Note that in the case of a background DOF we replace ψ h with 1 − ψ h . As E ( x i ) → x i become very small. If E = 0 then the contribution to the PUFEMsolution is null and the system matrix becomes singular. Thus, we define the four setsof constrained nodes as follows: X ve = { x i ∈ N ve | x i ∈ Γ h ff } , (12a) X pe = { x i ∈ N pe | x i ∈ Γ h ff } , (12b) X vb = { x i ∈ N vb | x i ∈ ¯Ω hcut , x i / ∈ Γ h ff } , (12c) X pb = { x i ∈ N pb | E ( x i ) < (cid:15) } . (12d)In the case of the constraints on Γ ff , i.e. (12a) and (12b), we have not observedparticularly small values of the E metric. However, we decided to include these nodes6n order to obtain a smoother solution for the embedded mesh. ¯Ω hcut refers to the areaof the background domain comprised of elements completely covered by the embeddedmesh. The definition in (12c) is meant to prevent any ill-conditioning, due to havinga non-unique solution. The second condition, that the background nodes do not lieon fluid-fluid interface, is meant to prevent a circular definition (e.g. v b ( x i ) = v e ( x i ),therefore v e ( x i ) = v b ( x i )). In the case of (12d), we chose to reduce the number ofconstrained DOF compared to X vb and generally assigned (cid:15) to be 0.1. In tests whichare not shown here, we saw that increasing the number of constrained pressure nodescan lead to a deterioration of the incompressibility condition at the domain level. Whilethis effect can be lowered by decreasing (cid:15) , the trade off is an increase in the system’scondition number. Following the definition in Eq. (8), (10) and (9), the discrete fields v h and p h can beexpanded into the sum of weighted basis functions: v h ( x ) = (cid:88) I (cid:2) − ψ h ( x ) (cid:3) ˆ φ Iv ( x )˜ v Ib + (cid:88) J ψ h ( x )ˆ ζ Jv ( x )˜ v Je , (13a) p h ( x ) = (cid:88) K (cid:2) − ψ h ( x ) (cid:3) ˆ φ Kp ( x )˜ p Kb + (cid:88) L ψ h ( x )ˆ ζ Lp ( x )˜ p Le . (13b)ˆ φ and ˆ ζ denote the piecewise polynomial functions defined on global and embeddedmeshes, respectively. Their subscripts are used in order to differentiate between first order( p ) and second order ( v ) functions, while the superscript marks the degree of freedom(DOF) index. Thus, the four sets of basis functions can be written as: { ˆ φ Iv } ≤ I ≤ N v,b , { ˆ φ Kp } ≤ K ≤ N p,b , { ˆ ζ Jv } ≤ J ≤ N v,e and { ˆ ζ Lp } ≤ L ≤ N p,e . N v,b , N v,e , N p,b and N p,e are the totalnumbers of DOFs corresponding to the V hb , V he , W hb and W he spaces. The nodal DOFvectors ˜ v Ib and ˜ v Je contain d entries each (e.g. ˜ v Ib = (cid:104) v Ib, , · · · , v Ib,d (cid:105) T ). Note that, bydesign, this definition varies according to where we evaluate the field. Thus, for x ∈ Ω he we have the definition in Equation (13a), while for x ∈ Ω hb \ Ω he the expansion simplifiesto: v h ( x ) = (cid:88) I ˆ φ Iv ( x )˜ v Ib and p h ( x ) = (cid:88) K ˆ φ Kp ( x )˜ p Kb . We can re-arrange the DOFs into four distinct vectors of unknown. For the back-ground mesh, for example, we have:˜ V b = (cid:104) v b, , . . . , v N v,b b, , . . . , v b,d , . . . , v N v,b b,d (cid:105) T and ˜ P b = (cid:104) p b , . . . p N p,b b (cid:105) T The ˜ V e and ˜ P e vectors are defined analogously.Algebraically, we can re-write the PUFEM weak form in Eq. (7) as: A bb A be B bb B be A eb A ee B eb B ee − B Tbb − B Teb − B Tbe − B Tbb ˜ V b ˜ V e ˜ P b ˜ P e = d b d e (14)7ere the A and B blocks result from the Laplacian and incompressibility terms, respec-tively. Their general structure takes the form: A αβ = A αβ . . .
00 0 A αβ and B αβ = B αβ ... B dαβ where α and β are placeholders for any b and e permutation. The inner sub-blocks aredefined as: ( A bb ) mn = (cid:90) Ω hb µ ∇ (cid:104) (1 − ψ h ) ˆ φ mv (cid:105) · ∇ (cid:104) (1 − ψ h ) ˆ φ nv (cid:105) d Ω (15a)( A be ) mn = (cid:90) Ω he µ ∇ (cid:104) (1 − ψ h ) ˆ φ mv (cid:105) · ∇ (cid:104) ψ h ˆ ζ nv (cid:105) d Ω (15b)( A eb ) mn = (cid:90) Ω he µ ∇ (cid:104) ψ h ˆ ζ mv (cid:105) · ∇ (cid:104) (1 − ψ h ) ˆ φ nv (cid:105) d Ω (15c)( A ee ) mn = (cid:90) Ω he µ ∇ (cid:104) ψ h ˆ ζ mv (cid:105) · ∇ (cid:104) ψ h ˆ ζ nv (cid:105) d Ω (15d)( B ibb ) mn = − (cid:90) Ω hb (cid:104) (1 − ψ h ) ˆ φ np (cid:105) ∂∂x i (cid:104) (1 − ψ h ) ˆ φ mv (cid:105) d Ω (15e)( B ibe ) mn = − (cid:90) Ω he (cid:104) ψ h ˆ ζ np (cid:105) ∂∂x i (cid:104) (1 − ψ h ) ˆ φ mv (cid:105) d Ω (15f)( B ieb ) mn = − (cid:90) Ω he (cid:104) (1 − ψ h ) ˆ φ np (cid:105) ∂∂x i (cid:104) ψ h ˆ ζ mv (cid:105) d Ω (15g)( B iee ) mn = − (cid:90) Ω he (cid:104) ψ h ˆ ζ np (cid:105) ∂∂x i (cid:104) ψ h ˆ ζ mv (cid:105) d Ω (15h)Vectors d g and d e are the right hand side vectors resulting from imposed Dirichletconditions. In order to apply the constraints on the system matrix, we need to distinguishbetween the Dirichlet condition and the constraints defined in Section 2.2.1. While theformer is trivial, in the case of the latter, the original system matrix line is replaced witha set of local basis functions evaluated at the constrained coordinate. For example, if wewanted to constrain an arbitrary background velocity degree of freedom, v mb,n , with thespatial coordinate x b,n , the new line of the system of equation would be equivalent to: v mb,n − (cid:88) i ∈ X m ˆ ζ iv ( x b,n ) v ib,n = 0 (16)All other nodes, whether background or embedded, velocity or pressure, are fixed analo-gously.We use X m to denote the set of nodes of embedded element in which we find x b,n .We use X m to denote the set of nodes of embedded element in which we find x b,n . Inthe next section, we expound the process involve in computing the weighted weak formintegral defined in (15). We will start by illustrating this process with A bb . The definition of this block, asdefined in Eq. (15a) can be rewritten as the sum of a classic weak form matrix, defined8n the background mesh, and a mixed classic-PUFEM matrix defined on the overlap:( A bb ) mn = (cid:90) Ω hb µ ∇ ˆ φ mv : ∇ ˆ φ nv d Ω+ (cid:90) Ω he − µ ∇ ˆ φ mv : ∇ ˆ φ nv + µ ∇ (cid:104)(cid:0) − ψ h (cid:1) ˆ φ mv (cid:105) : ∇ (cid:104)(cid:0) − ψ h (cid:1) ˆ φ nv (cid:105) (17)Thus we can separate the construction of A bb into two stages: first building a Laplacianblock using classic FEM and second subtracting the classic weak form from the overlaparea and replacing it with the PUFEM version. While the former is trivial, the later ismore challenging as the fields and integration area are defined on non-matching meshes.In order to overcome this, we construct an intersection or tertiary mesh of sub-elementswhich allows us relate the background and embedded meshes. Thus, at the start ofthe second stage, for each background element we want to identify all potentially inter-secting embedded elements. This can be achieved optimally by adapting one of severalapproaches found in the literature [59, 60, 61]. We then test the pairs of background andembedded candidates for intersection and we identify the area of overlap, see Fig. 3. Ifthe area is non-zero, then the resulting intersection will be a convex polygon with up to6 sides. Each area is subdivided into up to 4 subelements using a Delauney triangulationalgorithm and these are stored in the tertiary mesh. The quality of this new mesh is notrelevant since it is used for integration exclusively. To finalize the process, we loop overthe intersection mesh and we update the block by subtracting the overlap classic weakform and replacing it with the PUFEM one.Both A be and A eb are defined exclusively on Ω e and thus do not require the firststage introduced A bb . Since both are defined using fields from both meshes, the blocksare constructed by looping over the intersection mesh. In the case of A ee , the block canbuilt using the classic approach due to matching topologies for the space and fields.The B blocks are generated following an analogous procedure.
3. A PUFEM approach to ALE Navier-Stokes and FSI problems
Building on the foundations of Section 2, in this section, we present the applicationof PUFEM in the case of more complex settings, namely for ALE Navier-Stokes andFSI problems. Similarly to the Stokes flow case, we begin with a short review of theclassic ALE approach (Section 3.1), followed by a presentation of the PUFEM version(Section 3.3), where the focus is on how the method is adapted to deal with having atransient overlap between the two meshes. Finally, in Section 3.5 we outline a methodof coupling the fluid solver with a hyper-elastic solid problem.
Let us consider a time dependent physical domain Ω f ⊂ R d × [0 , T ]. At each time point t ∈ [0 , T ], we use the Ω f ( t ) notation to refer to the current spatial configuration of thefluid domain and let boundary Γ( t ) be its boundary. We also define the reference domainΛ f ⊂ R d which can be bijectively mapped to Ω f ( t ), for all t ∈ [0 , T ]. Hence, P : Λ f × [0 , T ] → Ω f denotes a mapping function used to relate the reference and spatial domains.9he non-conservative incompressible Navier-Stokes equation in arbitrary Lagrangian-Eulerian (ALE) form can be written as follows [11, 62]: find ( v , p ) such that ρ∂ t v + ρ ( v − ˆ v ) · ∇ v − ∇ · σ f = in Ω f ( t ) , (18a) ∇ · v = 0 in Ω f ( t ) , (18b) v = v D on Γ Df ( t ) , (18c) v = v S on Γ Dfs ( t ) , (18d) σ f · n = t N on Γ Nf ( t ) , (18e) v ( · ,
0) = v in Ω f (0) , (18f)for all t ∈ [0 , T ]. Compared to the Stokes flow in (1), the fluid model is augmented witha series of additional elements. The stress tensor is replaced by the symmetric Cauchystress tensor: σ f ( v ) = µ (cid:0) ∇ v + ∇ v T (cid:1) − p I (19)Inertial effects are introduced using the time derivative and nonlinear advection terms.The ∂ t ( · ) operator is the time derivative with respect to a fixed point in Λ f . The ALEadvection term is used to account for the arbitrary motion of the domain, where thearbitrary domain velocity field is defined as ˆ v = ∂ t P [62, 11]. Finally, we continue toassume that t n = .Using BDF(2) discretization, the classic non-conservative FEM problem for equa-tion (18) at a given time step n can be written as [62]: find ( v hn , p hn ) ∈ V hD × W h suchthat ∀ ( w h , q h ) ∈ V h × W h we have (cid:90) Ω hf ( t n ) ρ (cid:20) v hn − v hn − + v hn − t + (cid:16) v hn − ˆ v hn (cid:17) · ∇ v hn (cid:21) · w h d Ω+ (cid:90) Ω hf ( t n ) σ f ( v h ) : ∇ w h + q h ∇ · v hn d Ω =0 (20)Here, we consider the v hn − , v hn − ∈ V h to be the velocity fields at times t n − and t n − ,respectively, which have been transported via ALE mapping into the current spatialconfiguration.Due to the non-linearity of the system, we approximate the solution using a Newton-Raphson (NR) based algorithm [63, 64]. The problem is solved semi-monolithically inthat the mesh velocity field, ˆ v hn , is solved separately. If v s is provided analytically oris obtained from a weakly-coupled solid solver, then this process can be done at thebeginning of each time step. If the fluid and solid solvers are coupled in a monolithicsystem, then the mesh velocity is recomputed prior to each NR iteration in order toaccount for the current estimate of the solid deformation. The deformation of the mesh can be determined by an arbitrary problem which isguided by the known deformations on portions of the boundary. It is also well knownthat the choice of the problem can significantly impact mesh quality over time. Forthis reason, we propose to propagate the deformation of the boundary onto the entire10esh by solving a solid mechanics problem based on the nearly-incompressible neo-Hookean material model. In order to prevent an excessive deterioration of the mesh weconsider the material to be heterogeneous as we allow each element to stiffen in a squaredinverse relationship with a metric of its relative element quality. A similar idea has beensuccessfully implemented in [65].Let us consider the case where we want to compute the mesh deformation for time n . We choose Λ f to be the reference undeformed domain. Thus the solid problem canbe written as follows: find the arbitrary domain deformation field ˆ u such that ∇ · P g = in Λ f , (21a)ˆ u = ˆ u D on Γ Dref , (21b) P g · N = on Γ Nref . (21c)This relates to the domain deformation velocity as ˆ v = ∂ t ˆ u . Here, the (0) subscriptindicates that the divergence operator is defined in the Lagrangian coordinates and thefirst Piola-Kirchhoff tensor ( P g ) is given as: P g = µ g (cid:20) J (cid:18) F − F : F F − T (cid:19) + κJ ( J − F − T (cid:21) (22)where F = ∇ ˆ u is the deformation gradient and J = det( F ). The heterogeneous stiffnessparameter , µ g , is defined as: µ g ( τ ) = µ ref (cid:18) Q ( τ ) Q ( τ ) (cid:19) where Q ( τ ) = rR Here µ ref is the reference stiffness; Q : T → R + is an element quality metric defined asthe ratio between the incircle ( r ) and circumcircle ( R ) of input triangular element; τ and τ denote the original and current shape of the element, respectively. As the qualityof the element deteriorates, Q ( τ ) decreases, making it stiffer and, thus, less likely tofurther deform. κ = 10 is a penalty factor on volumetric deformation. This choice wasdone heuristically and it aims to penalize excessive decreases in element size.In the discrete time setting, we chose to use the previous time step configuration as thereference, undeformed domain. For the normalization factor Q ( τ ), we use the elementat time zero as our reference. Given a discrete time point t n , we use the deformationfield to compute the current space mapping function follows: P h ( X h , t n ) = P h ( X h , t n − ) + ˆ u hn ( X h ) (23)Based on this, the mesh velocity is defined as ˆ v hn = ˆ u hn ∆ t . In Section 2.2 we discussed the setup of PUFEM for flow around a rigid solid. Nowwe expand these original concepts in order to treat the case where the two meshes areallowed to move independently from each other.A key requirement for the extension of PUFEM to ALE is the treatment of thematerial time derivative. Let Λ hb and Λ he be the two reference domains for the background11nd embedded components. P b : Λ hb × [0 , T ] → Ω hb and P e : Λ he × [0 , T ] → Ω he are thetwo arbitrary mapping functions, where Ω hb , Ω he ∈ R d × [0 , T ] are the respective physicaldomains. We also introduce the ∂ bt ( · ) and ∂ et ( · ) as the partial time derivatives withrespect to fixed points the global and embedded reference domains, respectively. Forease, we consider the ALE expansion of the material time derivative on an arbitrarysemi-discrete (in space) scalar field f h ( t ) ∈ V h which can be written in PUFEM formas f h = (1 − ψ h ) f hb + ψ h f he , where f hk ( t ) ∈ V hk, ∗ , k ∈ { e, b } . Based on this, we can splitthe time material derivative of f into three terms: Df h Dt = D (cid:2) (1 − ψ h ) f hb + ψ h f he (cid:3) Dt = (1 − ψ h ) Df hb Dt + ψ h Df he Dt + ( f he − f hb ) Dψ h Dt Furthermore, we choose Λ hb to be the reference domain of f hb and Λ he for f he and ψ h .Thus, we can continue to expand the first term into its ALE form such that: Df hb Dt (cid:12)(cid:12)(cid:12)(cid:12) x = ∂f hb ∂t (cid:12)(cid:12)(cid:12)(cid:12) x + v h · ∇ f hb = (cid:18) ∂ bt f hb − ∂ P hb ∂t · ∇ f hb (cid:19) + v h · ∇ f hb = ∂ bt f hb + ( v h − ˆ v hb ) · ∇ f hb Similarly, we also get: Df he Dt (cid:12)(cid:12)(cid:12)(cid:12) x = ∂ et f he + (cid:16) v h − ˆ v he (cid:17) · ∇ f hb and Dψ h Dt (cid:12)(cid:12)(cid:12)(cid:12) x = (cid:16) v h − ˆ v he (cid:17) · ∇ ψ h In the last case, we assumed that ∂ et ψ h = 0. The total material time derivative can thusbe expressed as: Df h Dt = (1 − ψ h ) ∂ bt f hb + ψ h ∂ et f he + v h · ∇ f h − (1 − ψ h )ˆ v hb · ∇ f hb − ψ h ˆ v he · ∇ f he − ( f he − f hb )(ˆ v he · ∇ ψ h )Moving to the discrete form, let us define Ω hb,n and Ω he,n , the global and embeddedmeshes, respectively, for a given time step n . Given that X vb , X ve , X pb and X pe changewith the alignment of the meshes with respect to each other at each time step, let usdefine V hD,n , V h ,n and W hn as the test and trial spaces at discrete time t n . Thus thenew discrete weak form of the ALE Navier-Stokes problem can be written as follows:find ( v hn , p hn ) ∈ V hD,n × W hn such that for any ( w h , q h ) ∈ V h ,n × W hn we have that (cid:90) Ω hb ( t n ) ρ (cid:18) v hn − v hn − + v hn − t + v hn · ∇ v hn (cid:19) · w h d Ω+ (cid:90) Ω hb ( t n ) − ρ (cid:104) (1 − ψ hn )(ˆ v hb,n · ∇ ) v hb,n + ψ hn (ˆ v he,n · ∇ ) v he (cid:105) · w h d Ω+ (cid:90) Ω hb ( t n ) − ρ (ˆ v he · ∇ ψ hn )( v he,n − v hb,n ) · w h + σ hf,n : ∇ w h − q h ∇ · v hn d Ω =0 (24)where v hn = (1 − ψ h ) v hb,n + ψ h v he,n . Both partial derivatives ∂ gt ( · ) and ∂ et ( · ) have beenapproximated using the BDF(2) scheme. In order to solve this system of equations, we12gain employ the technique described in [64] and [63]. The structure of the resultingJacobian matrix is largely based on that of system matrix described in Eq. (14). Thechanges arise in the sub-block A where in addition to the Laplacian we also have massand advection terms resulting from linearisation of the residual: (cid:16) A ijbb (cid:17) MN = (cid:90) Ω hb,n ρδ ij (cid:20) t ψ hb ψ hb ˆ φ Mv ˆ φ Nv + ψ hb ˆ φ Mv (cid:0) v hn · ∇ (cid:1) (cid:16) ψ hb ˆ φ Nv (cid:17)(cid:21) d Ω+ (cid:90) Ω hb,n ρψ hb ψ hb (cid:34) ˆ φ Mv ˆ φ Nv ∂ (cid:0) v hn (cid:1) i ∂x j − δ ij ˆ φ Mv (cid:16) ˆ v hb,n · ∇ (cid:17) ˆ φ hv (cid:35) d Ω − (cid:90) Ω hb,n ρδ ij ˆ φ Mv ˆ φ Nv (cid:16) ˆ v he · ∇ ψ hb (cid:17) − µδ ij ∇ (cid:16) ψ hb ˆ φ Mv (cid:17) : ∇ (cid:16) ψ hb ˆ φ Nv (cid:17) d Ω+ (cid:90) Ω hb,n µ ∂ψ hb ˆ φ Mv ∂x j ∂ψ hb ˆ φ Nv ∂x i d Ω (25)The other blocks resulting from permutations of background/embedded test and trialfunctions are obtained analogously. Here, i and j correspond to the spatial orientationof test and trial DOF, respectively. M and N represent the indexes of the mesh nodes.For ease, ψ hb is used as shorthand for (1 − ψ h ). The application of Dirichlet boundaryconditions and additional constraints is largely the same as the one described in Sec-tion 2.2. The exception is that while the Jacobian is re-used for multiple time iterations,the constraints such as v b = x e for x ∈ X vb need to be re-adjusted based on the newoverlap configuration. Moving the background and embedded meshes independently from each other intro-duces new challenges: having a different area of overlap in consecutive time steps and, byextension, having transient sets X vb and X pb . While the latter is of less importance, theformer can be significant factor of instability due to the need to use previous time stepsolutions in order to approximate the velocity time derivative, as shown in Eq. 24. Thus,if a node is fixed at one time t n and subsequently active at time t n +2 , our estimation ofthe acceleration is dependent on what we chose to assign t n . While, this concern has beenpartly addressed in Section 2.2.1, it only covers the case of background nodes which arefound under the fluid area of the embedded mesh. Thus, in this case we can assign thema meaningful value through interpolation. This leaves, however, the case of backgroundnodes found under the solid. Here there are two scenarios. In the first case our area ofenrichment is sufficiently thick and the time step is sufficiently small such that when abackground node, after it leaves solid, remains in the overlap area for at least two moretime steps. This would allow for sufficient time for the node to obtain meaningful data inorder to compute the derivative. However, this would be difficult to control, particularlyin FSI. An alternative solution that we propose is to constrain the value of backgroundnodes in the solid to be equal to the interpolated solid velocity field, which is generallycontinuous to that of the fluid flow. If this field is provided or is computed as part ofthe FSI solution, then this new task is trivial. If however, only the surface motion isprovided, one can potentially create a field by solving an artificial problem, such as adiffusion one. It should be noted though that this last approach is only theoretical andwas never implemented for this work. 13 .5. PUFEM for FSI To illustrate the efficacy of the PUFEM approach for FSI, in this section we elab-orate on a monolithic approach to coupling the ALE Navier-Stokes solver described inSection 3.3 with a generic quasi-static non-linear solid. Let Ω s, and Ω s,t ⊂ R d be ourreference and deformed (at time t ) solid domains. Thus, the strong form of the FSI andALE mesh problems can be written as: find ( v f , p f , u s , p s , ˆ u e ) such that ρ∂ t v f + ρ ( v f − ˆ v ) · ∇ v f − ∇ · σ f ∇ · v f v f σ f · nv f ( · ,
0) = = 0= v f,D = t f,N = v in Ω f , in Ω f , on Γ Df on Γ Nf , in Ω f, , (26) ∇ · P s ( u s , p s ) J − u s P s ( u s , p s ) · N = = 0= u s,D = t s,N in Ω s, in Ω s, on Γ Ds, on Γ Ns, (27) σ f · n + J − P s ( u s , p s ) · N v f = = ˙ u s on Γ fs on Γ fs (28) ∇ · P g (ˆ u e )ˆ u e P g (ˆ u e ) · N = = u f = in Ω e, on Γ fs , on Γ ff , (29)The FSI problem can be divided into four parts: the fluid problem in (26), the quasi-static solid problem in (27), the coupling conditions in (28) and the ALE mesh problemin (29). For illustrative purposes, we chose to use the neo-Hookean model for the solid,resulting in the following form of the first Piola-Kirchhoff stress tensor: P s = µ s J s (cid:20) F s − F s : F s F − Ts (cid:21) − p s J s F − Ts (30)where F s = ∇ u s and J s = det ( F s ). Note, we chose to simplify the solid problem byconsidering the loading process to be quasi-static. The constitutive law used for thearbitrary mesh deformation problem has been previously described in Section 3.2. ThePUFEM discrete setting is elaborated on in the next section. The main difference compared to the ALE Navier-Stokes set-up described in Sec-tion 3.3 is that now we also need to solve the deformation/translation of the solid, ratherthan it being given. In order to take advantage of the definition of the PUFEM spacesin (7), (10) and (9) which extend into both Ω hf and Ω hs , for each time step we computethe solid’s velocity v hn , rather than its displacement. The deformation is obtained usingthe backward Euler scheme ( u hs,n = u hs,n − + v hs,n ∆ t ). Thus, the problem coupling isensured by the fact that the fluid and solid velocity DOF are identical on Γ fs .14he fully discrete weak form can be written as: find ( v hn , p hn ) ∈ V hD,n × W hn such that ∀ ( w h , q h ) ∈ V h ,n × W hn (cid:90) Ω hf,n ρ (cid:18) v hn − v hn − + v hn − t + v hn · ∇ v hn (cid:19) · w h d Ω+ (cid:90) Ω hf,n − ρ (cid:104)(cid:0) − ψ hn (cid:1) ˆ v hb,n · ∇ v hb,n + ψ hn ˆ v he,n · ∇ v he,n (cid:105) · w h d Ω+ (cid:90) Ω hf,n − ρ (ˆ v he · ∇ ψ hn )( v he,n − v hb,n ) · w h d Ω+ (cid:90) Ω hf,n σ hf,n : ∇ w hn + q h ∇ · v hn d Ω+ (cid:90) Ω hs, P s : ∇ w h + q h ( J s − d Ω = 0 (31)The structure of the system of equations which has to be solved during each Newton-Raphson remains largely the same. Note, the pressure is allowed to be discontinuousacross Γ fs In order to examine the performance of PUFEM against the classic mixed FEMtechnique, we implemented both approaches in the MATLAB 2017b language. Themeshes used for the test were produced using the MESH2D package [66, 67]. All thesimulation where run on MacBook pro 2017 running on macOS Sierra, with an IntelCore i7-7820HQ processors and 16 GB of 2133 MHz LPDDR3 RAM.
4. Numerical Results
In order to illustrate the comparable accuracy of PUFEM to that of standard bound-ary fitted FEM, we consider a basic steady-state incompressible Stokes flow problem,where we compute the errors in the velocity and pressure fields on series of meshes.More specifically, we consider a total domain Ω f ∪ Ω s = [0 , where a rigid cylindricalobstacle (with a radius of 0.15) is placed in the centre, see Fig. 4. The inflow (left side) isconstrained using the Dirichlet condition, v inflow = [1; 0] T , which is constant throughoutthe boundary patch. On the side walls, we apply a reflection Dirichlet condition on thevelocity field, such that v y = 0, and on the outflow we impose a zero traction condition.Fluid viscosity was set to 1.0. The purpose of this arrangement is to concentrate theboundary layer flow to the vicinity of cylinder, which coincides with the area enrichedby the embedded mesh.To ease comparison, we try to use meshes with small variations in element size. Forboth approaches, we consider five refinement levels, with h ranging from approximately0 . . h g ≈ h e . The latter mesh is composed of two regions: the cylinder, whichis only used to define the solid surface, and the fluid enrichment area, which forms a15ing with a thickness 0.161. We shall refer to this grid as M1. More details about theirstatistics and that of the classic mesh set can be found in Table 1.The resulting velocity and pressure fields and their subcomponents can be seen inFig. 5. The PUFEM weighted sum result and the classic solution appear to be verysimilar. Additionally, there are no obvious non-physical behaviours, such as field jumps,which might have been generated by the coupling.The results of the convergence test are shown in Fig. 6. Note, that we did not relyon an analytical reference solution to compute the error. In order to obtain a reference,the flow field was computed using the classic approach and a very fine, boundary fittedgrid with an element size of h = 0 . a priori estimate in Eq. 6. We observe that the velocityproduces slightly suboptimal results, however, this is likely due to both projection errorsand non-conformity between solutions. We also notice that the difference in errors acrossthe approaches are very small. Given the mesh size similarity among the methods,PUFEM was not expected to exceed the performance of the classical approach. Thus, theaccuracy of the new method is ultimately capped by the local approximation capabilityof the background and embedded grids. Furthermore, the coupling, representing themain potential error source, appears to have little deteriorating impact. To examine PUFEM performance for steady-state Navier-Stokes, we consider theproblem outlined in Fig 4 for Reynolds numbers of 30 and 100. We keep using the [0 , domain with R = 0 .
15 cylindrical obstacle in the middle. The fluid parameters are setto: ρ = 1 . µ = 0 .
01. The boundary conditions are kept largely the same withthe exception of the inflow velocity. Thus, the x direction component is set to 1 . /
3, depending on whether we want to compute the Re
30 or 100 cases. Again, ourerrors estimates are based on flow results computed using the classic approach and a finegrid with h = 0 . h b ≈ h e . See Table 1 forthe corresponding mesh statistics.A set of sample results for the velocity and pressure at Re
100 can be seen in Fig. 7.Despite the increase in flow complexity, the quality of the PUFEM solution remainscomparable to that of the classic result.The error plots for the different cases are found in Fig. 8. A first observation is thatthe convergence behaviour for Re
30 and 100 remains relatively unchanged from thatobserved in Stokes flow problem. Furthermore, the differences in accuracy between clas-sic and PUFEM (M1) approaches remain small. As previously noted in the Stokes case,when using similar spatial discretization, we do not expect the new approach present asignificant advantage. It continues to perform well in these circumstances and is seem-ingly unaffected by the rise in Re . Switching to M2, we see an overall decrease in error.This suggest that PUFEM can potentially leverage local solution enrichment to improveoverall accuracy. While a similar behaviour can be replicated in the case of boundary16tted grids using adaptive meshing, the main strength of the new method resides in itsability to better withstand large deformations. In order to evaluate the stability of the PUFEM method in the case of time depen-dent problems and also to compare it to other published flow results, we considered theSch¨afer-Turek benchmark [68], case 2D-2 with Re = 100 (see Fig. 9).The inflow constraint is imposed on the velocity using the following Dirichlet condi-tion: v = (cid:20) Uy ( W − y ) W min( t, (cid:21) After linearly ramping up, the inflow profile remains constant and we run the problem for10 s in simulation time, enough to initiate the characteristic shedding and also to reachsteady state oscillations. Here U = 1 . m/s is the reference flow speed and W = 0 . m is the width of the tube. We also impose a no slip condition on the side walls and zerotraction on the outflow. In order to achieve the correct Reynolds number, density ( ρ )and viscosity ( µ ) are set to 1 kg/m and 0 . P a · s , respectively.In total, six solutions are produced using PUFEM based on two mesh sets (PU1 andPU2) and three time steps (i.e. 0.01, 0.002 and 0.001 s ). Each mesh set consists ofa regular background grid and a boundary fitted embedded grid. The fluid region ofthe latter forms a ring with a thickness of 0.1 m around the cylinder, see Fig. 9. Forsimplicity, each set uses a constant h for both background and embedded grids, with thevalue set to 0.025 and 0.0125 cm in the case of PU1 and PU2, respectively. Additionalmesh statistics are found in Table 2. To aid the comparison, we have also produced aset of flow results using the classical approach using a quasi-regular grid with h = 0 . s .Fig. 10 shows an example set of flow results (i.e. velocity components and pressure)computed using classic and PUFEM approaches. In the top part, we have a snapshotof the flow at one of the time points of maximum lift (i.e. 9.57 s in simulation timefor PUFEM). The results are split into subcomponents corresponding to the boundaryfitted mesh, background, embedded meshes and the weighted sum. In the second halfof the figure, we include six velocity magnitude snapshots at the different stages of theshedding process. Each frame is split in two, with the right side showing the PUFEMfield and the left displaying the classic, showing the wake patterns are in good agreement.Furthermore, the transition between the embedded and background grids retains itssmoothness, despite the transient behaviour of the solution.From a quantitative perspective, we also computed the coefficients of drag ( c d ), lift( c l ), Strouhal number( St ) and pressure drop across the cylinder (∆ p ), defined as: c d = 2 F d ρ ¯ U D , F d = (cid:90) Γ fs µ ∂v t ∂n n y − pn x d Γ c l = 2 F l ρ ¯ U D , F l = − (cid:90) Γ fs µ ∂v t ∂n n x + pn y d Γ St = f D ¯ U , D is the diameter and v t is the tangential fluid velocity. These values have beencompiled in Table 3, which contains all the PUFEM and classic results, as well as thosefrom literature. For similar space and time discretization, i.e. entries two and three ofthe table, the new and classic approaches provide similar estimations of the with errorsof: 0 .
03% for c d , 0 .
22% for c l , 3% for St , and 0 .
16% for ∆ p . By comparing to theliterature values [68], we see that the PUFEM estimations of St and ∆ p are converginginside the recommended interval. However, c d and c l are generally underestimated, withminimum errors of 1 .
18% and 1 . s were approx. 1.9 and 2.9 hours, respectively. In this case, where the mesh overlapis fixed, the main factor responsible PUFEM’s lag appears to be the residual’s assembly,with an average running time 0.42 s , compared to 0.097 s for the classic. While notaddressed in this work, there are multiple avenues which could be used to reduce the thisdisparity: restricting the interface mesh to the area where ψ is strictly between 0 and 1,contracting this area to a thinner strip close to Γ h ff and optimizing the quadrature rulebeing used, to name a few. Moving towards moving domain applications, we define a simple ALE benchmarkproblem. Using a similar to setup to that of the Turek benchmark, we now place thecylinder in the middle of the domain, and we oscillate it along the long axis of the tube,while we allow for free outflow at the left and right ends. The cylinder displacement asa function of time is defined as: d ( t ) = (cid:26) A sin( ωt ) ω − At cos( ωt ) , t ≤ − A cos( ωt ) , t > A is the maximum amplitude and ω is the frequency. We chose to set A = 0 . ω = 2 π s − , resulting in a maximum velocity of approx 1 . ms − and a Re ≈ cut .For the spatial discretization, we define a new ALE-PU1 mesh set to be consistentwith the new geometry (i.e. in accordance to the repositioned cylinder) and with identicalmesh statistics as PU1 (i.e. due to recycled connectivity matrices). For the classical ap-proach, a new boundary fitted mesh (ALE-FEM1) was constructed with 14708 elementsand h = 2 . t = 0 . s was used for the temporal discretization.For a qualitative comparison between methods, we maintained the same format fordisplaying the flow results as in the case of the Sch¨afer-Turek benchmark, see Fig. 11.The top part of the figure is a snapshot at t = 2 s showing the velocity componentsand pressure field as computed on the boundary fitted, the background and embedded18rids. The corresponding PUFEM total field is also included. The middle part displayssnapshots of the velocity magnitude at different stages of a transition from right to left.In general, the two methods seem to be in agreement, with only minor differences. Incontrast to previous examples, the ALE problem is the first case in this work whichincludes the concept of transient fixed nodes, see Fig. 12 .Thus, the good quality ofour solution now also suggests the effectiveness of using the interpolation strategy onproviding nodes with a meaningful solution while they are temporarily fixed.We also computed the mean surface drag force ( F d ) and pressure across the left andright ends of the cylinder (∆ p ). Fig. 11 (bottom) shows the values of these coefficients atall time points of the simulation. The maximum overall errors that we obtain are 4 .
0% forthe mean F d and 0 .
8% for ∆ p . Potential sources which may influence the error (excludingPUFEM) include mesh deterioration in the case of the classic approach and also the factthat we did not aim to minimize the errors through spatial and temporal refinement.Despite this, these results suggest that PUFEM can be used in an ALE context toobtain reasonably accurate estimations of surface stresses and pressure experienced bythe solid. This quality is particularly relevant in the case FSI, where the balance oftractions is part of the coupling conditions. To examine the suitability of PUFEM for FSI applications, we present an example ofan idealised 2D aortic valve, see Fig. 13. The FSI domain consists of two components:one half of the lumen and one leaflet. In order to reduce computation cost, we assumeto have rigid aortic walls and symmetry along the long axis. We also do not includecoronary flow. The fluid is modelled as an incompressible Navier-Stokes fluid, while thesolid is assumed to behave as an incompressible neo-Hookean material [69]. The maindomain characteristics and constitutive law parameters can be found in Table 4. At theinflow we impose a single pulse with parabolic a profile: v ( y, t ) = (cid:20) V max y (2 W − y ) W sin ( πt )0 (cid:21) (33)where V max = 0 . m/s is the maximum velocity, W is the width of the inflow and y is thespatial coordinate. The total in-simulation duration is T = 1 s , or one pulse, which wediscretise using constant time steps of ∆ t = 0 . s . Thus we can observe, the flow fieldsboth during valve deflection as well as relaxation. On the surface of the wall and leafletwe impose a no-slip condition, while on the long axis we apply a symmetry condition.For the outflow, we impose a zero traction Neumann condition.In order to reduce the element deterioration and the loss of accuracy which mayderive from it, the mesh used in the classic approach is built as follows. First, a mesh isconstructed on the reference domain and we run the simulation up to the t mid = 0 . s .The time point was chosen heuristically in order to capture the configuration of thedomain during mid valve deflection. This data is then used to construct a second mesh.Using the connectivity matrix of this grid, a final mesh is defined by interpolating the t = 0 . s node coordinates onto the space at t = 0 s . Thus, the classic grid is comprisedof 8248 elements, with mean and maximum h values of 0 . . mm , respectively.For PUFEM, the two grids are generated in the initial domain configuration. Weused a background mesh with 4380 element and an average h of 0.7 mm. The embeddedgrid is formed of 844 elements with mean h of 0.5 mm.19he main flow results are summarized in Fig. 14: the velocity magnitude fields atdifferent stages of the pulse (A); the grid deformation at half-way deflection and endtime point (B); and the element quality distribution as it evolves in time (C). From aqualitative perspective, the flow physics appear to be quite similar at all time stages.While significant differences appear in the second part of the simulation, we believe thatthese are generally due to a combination of the grid deterioration in the case of theclassic approach and decelerating flows, which are more sensitive to its effects. In fact,without the use midway meshing, to try to attenuate these effects, we observed muchmore striking difference between methods (results not shown here). The deteriorationof grid quality for the classic approach can be clearly seen in (B). Thus, while at T / T (end of relaxation) many of the elements close to the tip of the leaflet appear tobe squeezed against the axis of symmetry. On the other hand PUFEM does not showsuch clear signs, suggesting that mesh can undergo even more significant deformation.This is also reflected in (C), where for PUFEM’s embedded grid the minimum elementquality only temporarily drops bellow 0.5. In the case of the classic approach, the elementpopulation is selected from the interval [1 . , .
7] cm on the x axis, as shown in (B). Herewe see that while a majority number of elements retains an excellent quality, for manythe value is very low, close to collapsing. The reason for having a much wider populationspread at the start and end of the simulation is a result of using the midway meshingstrategy. In return, the population appears to be much more skewed towards one in themiddle half of the simulation, when the we experience the highest inflow velocity.From the solid mechanics perspective, a good agreement can also be observed for thevalve tip deflection, see Fig. 15, with maximum absolute errors of 0 .
21 and 0 . mm in x and y directions, respectively.Combined, these preliminary results on FSI suggest not only show a good level ofaccuracy for PUFEM, comparable to that of the classic approach, but also to a greaterdegree of flexibility, where we can better retain element quality while undergoing signif-icant mesh deformation.
5. Conclusion
In this work we have presented a novel approach to solving fluid flow and FSI problemsbased on the partition of unity (PUFEM). Thus, the total flow solution is defined asweighted sum of background and embedded fields. The weighting field, which is definedon the embedded mesh, is used to facilitate a smooth transition between areas where thesolution is dominated by one mesh or the other.The accuracy and stability of the method has been assessed on a range of 2D tests withan increasing level complexity. For comparison, all tests were replicated using a classicboundary fitted FEM approach which served as our golden standard. Using steady-state results for flow around the cylinder, we showed that PUFEM displayed nearlyoptimal convergence rates for both Stokes and Navier-Stokes, and that the the accuracyof the two approaches is nearly indistinguishable for similar spatial discretization. Thissuggests that PUFEM can leverage the approximation power of both meshes with littledeterioration as a result of the coupling process.In the case of the Sch¨afer-Turek, the flow results produced by PUFEM where shownto be in good agreement with the classic approach. Additionally, the estimation of surface20tresses, pressure and wake frequency appeared to be in accord with literature values.To assess the accuracy and stability of PUFEM in the ALE setting, we designeda problem where flow is generated in an open-ended channel by an oscillating, rigidcylinder. The method appeared robust to the changing of the overlapping area andthe set of constrained nodes. Additionally, a good agreement was maintained betweenapproaches in the estimation of surface drag forces and pressure drop.An idealised aortic valve FSI problem was chosen for testing of realistic applications.Given the set up the problem, the resulting flow is complex with the occurrence ofrecirculation and small vortexes during relaxation. In addition to displaying a high levelof similarity between the physical behaviour obtained using the two approaches, we whereable to gauge the evolution of mesh quality over the course of the simulation. Thus, theuse of overlapping meshes appears to produce the desired effect of reducing the wellknown mesh degradation associated with large deformations.Combined, these results indicate the potential applicability of PUFEM to a widerange of FSI applications, allowing for boundary fitted fluid-solid interface and targetedsolution space enrichment, in the absence of the need for re-meshing or the need foruser-defined stabilization parameters.Future work will focus on providing a theoretical analysis on the stability and conver-gence of the PUFEM solver. Other extensions of this work will include porting it into aparallel 3D implementation [70] and the study of different approaches to include contactmechanics.
Acknowledgments
D.N. acknowledges funding form the Engineering and Physical Sciences (EP/N011554/1and EP/R003866/1). A.M. gratefully acknowledges financial support from the SwedishResearch Council under Starting Grant 2017-05038 and from the Wenner-Gren foun-dation under travel grant SSh2017-0013. JH acknowledges the financial support of theSwedish Research Council under Grant 2018-04854. This work is funded by the King’sCollege London and Imperial College London EPSRC Centre for Doctoral Training inMedical Imaging (EP/L015226/1). This work is supported by the Wellcome EPSRCCentre for Medical Engineering at King’s College London (WT 203148/Z/16/Z) and bythe National Institute for Health Research (NIHR) Biomedical Research Centre award toGuy and St Thomas’ NHS Foundation Trust in partnership with King’s College London.The views expressed are those of the authors and not necessarily those of the NHS, theNIHR or the Department of Health.
Additional Information
Declarations of interest: none. 21 eferencesReferences [1] K. Stein, R. Benney, V. Kalro, A. Johnson, T. Tezduyar, K. Stein, R. Benney, V. Kalro, A. Johnson,T. Tezduyar, Parallel computation of parachute fluid-structure interactions, in: 14th AerodynamicDecelerator Systems Technology Conference, 1997, p. 1505.[2] V. Kalro, T. E. Tezduyar, A parallel 3D computational method for fluid–structure interactions inparachute systems, Computer Methods in Applied Mechanics and Engineering 190 (3-4) (2000)321–332.[3] K. Takizawa, T. Spielman, T. E. Tezduyar, Space–time FSI modeling and dynamical analysis ofspacecraft parachutes and parachute clusters, Computational Mechanics 48 (3) (2011) 345.[4] Y. Bazilevs, M.-C. Hsu, I. Akkerman, S. Wright, K. Takizawa, B. Henicke, T. Spielman, T. Tez-duyar, 3D simulation of wind turbine rotors at full scale. Part I: Geometry modelling and aerody-namics, International Journal for Numerical Methods in Fluids 65 (1-3) (2011) 207–235.[5] Y. Bazilevs, M.-C. Hsu, J. Kiendl, R. W¨uchner, K.-U. Bletzinger, 3D simulation of wind turbine ro-tors at full scale. Part II: Fluid–structure interaction modelling with composite blades, InternationalJournal for Numerical Methods in Fluids 65 (1-3) (2011) 236–253.[6] M.-C. Hsu, D. Kamensky, F. Xu, J. Kiendl, C. Wang, M. C. Wu, J. Mineroff, A. Reali, Y. Bazilevs,M. S. Sacks, Dynamic and fluid–structure interaction simulations of bioprosthetic heart valves usingparametric design with T-splines and Fung-type material models, Computational mechanics 55 (6)(2015) 1211–1225.[7] H. Gao, L. Feng, N. Qi, C. Berry, B. E. Griffith, X. Luo, A coupled mitral valve–left ventricle modelwith fluid–structure interaction, Medical engineering & physics 47 (2017) 128–136.[8] M. McCormick, D. Nordsletten, P. Lamata, N. P. Smith, Computational analysis of the importanceof flow synchrony for cardiac ventricular assist devices, Computers in biology and medicine 49(2014) 83–94.[9] K. Lau, V. Diaz, P. Scambler, G. Burriesci, Mitral valve dynamics in structural and fluid–structureinteraction models, Medical engineering & physics 32 (9) (2010) 1057–1064.[10] B. Baillargeon, I. Costa, J. R. Leach, L. C. Lee, M. Genet, A. Toutain, J. F. Wenk, M. K. Rausch,N. Rebelo, G. Acevedo-Bolton, et al., Human cardiac function simulator for the optimal design of anovel annuloplasty ring with a sub-valvular element for correction of ischemic mitral regurgitation,Cardiovascular engineering and technology 6 (2) (2015) 105–116.[11] Y. Bazilevs, K. Takizawa, T. E. Tezduyar, Computational fluid-structure interaction: methods andapplications, John Wiley & Sons, 2013.[12] A. Verkaik, M. Hulsen, A. Bogaerds, F. van de Vosse, An overlapping domain technique couplingspectral and finite elements for fluid flow, Computers & Fluids 100 (2014) 336–346.[13] C. W. Hirt, A. A. Amsden, J. Cook, An arbitrary Lagrangian-Eulerian computing method for allflow speeds, Journal of computational physics 14 (3) (1974) 227–253.[14] J. Donea, S. Giuliani, J.-P. Halleux, An arbitrary Lagrangian-Eulerian finite element method fortransient dynamic fluid-structure interactions, Computer methods in applied mechanics and engi-neering 33 (1-3) (1982) 689–723.[15] T. J. Hughes, W. K. Liu, T. K. Zimmermann, Lagrangian–Eulerian finite element formulation forincompressible viscous flows, Computer methods in applied mechanics and engineering 29 (3) (1981)329–349.[16] R. Van Loon, P. Anderson, F. Van de Vosse, S. Sherwin, Comparison of various fluid–structureinteraction methods for deformable bodies, Computers & structures 85 (11-14) (2007) 833–843.[17] T. E. Tezduyar, S. Sathe, Modelling of fluid–structure interactions with the space–time finite ele-ments: solution techniques, International Journal for Numerical Methods in Fluids 54 (6-8) (2007)855–900.[18] A. M. Bavo, G. Rocatello, F. Iannaccone, J. Degroote, J. Vierendeels, P. Segers, Fluid-structureinteraction simulation of prosthetic aortic valves: comparison between immersed boundary and arbi-trary Lagrangian-Eulerian techniques for the mesh representation, PloS one 11 (4) (2016) e0154517.[19] R. Glowinski, T.-W. Pan, J. Periaux, A fictitious domain method for Dirichlet problem and appli-cations, Computer Methods in Applied Mechanics and Engineering 111 (3-4) (1994) 283–303.[20] R. Glowinski, T.-W. Pan, T. I. Hesla, D. D. Joseph, A distributed Lagrange multiplier/fictitiousdomain method for particulate flows, International Journal of Multiphase Flow 25 (5) (1999) 755–794.
21] I. Babuˇska, The finite element method with Lagrangian multipliers, Numerische Mathematik 20 (3)(1973) 179–192.[22] L. Zhang, A. Gerstenberger, X. Wang, W. K. Liu, Immersed finite element method, ComputerMethods in Applied Mechanics and Engineering 193 (21-22) (2004) 2051–2067.[23] L. Zhang, M. Gay, Immersed finite element method for fluid-structure interactions, J. Fluids Struct.23 (6) (2007) 839–857.[24] D. Boffi, N. Cavallini, L. Gastaldi, The finite element immersed boundary method with distributedlagrange multiplier, SIAM J. Numer. Anal. 53 (6) (2015) 2584–2604.[25] D. Boffi, L. Gastaldi, A fictitious domain approach with lagrange multiplier for fluid-structureinteractions, Numer. Math. 135 (3) (2017) 711–732.[26] A. J. Gil, A. A. Carre˜no, J. Bonet, O. Hassan, The immersed structural potential method forhaemodynamic applications, Journal of Computational Physics 229 (22) (2010) 8613–8641.[27] N. D. Dos Santos, J.-F. Gerbeau, J.-F. Bourgat, A partitioned fluid–structure algorithm for elasticthin valves with contact, Computer Methods in Applied Mechanics and Engineering 197 (19-20)(2008) 1750–1761.[28] D. Kamensky, M.-C. Hsu, D. Schillinger, J. A. Evans, A. Aggarwal, Y. Bazilevs, M. S. Sacks, T. J.Hughes, An immersogeometric variational framework for fluid–structure interaction: Applicationto bioprosthetic heart valves, Computer methods in applied mechanics and engineering 284 (2015)1005–1053.[29] R. Van Loon, P. D. Anderson, J. De Hart, F. P. Baaijens, A combined fictitious domain/adaptivemeshing method for fluid–structure interaction in heart valves, International Journal for NumericalMethods in Fluids 46 (5) (2004) 533–544.[30] A. Gerstenberger, W. A. Wall, An extended finite element method/Lagrange multiplier based ap-proach for fluid–structure interaction, Computer Methods in Applied Mechanics and Engineering197 (19-20) (2008) 1699–1714.[31] F. Alauzet, B. Fabr`eges, M. A. Fern´andez, M. Landajuela, Nitsche-XFEM for the coupling of an in-compressible fluid with immersed thin-walled structures, Computer Methods in Applied Mechanicsand Engineering 301 (2016) 300–335.[32] M. Landajuela, Coupling schemes and unfitted mesh methods for fluid-structure interaction, Ph.D.thesis, Universit´e Pierre et Marie Curie, Paris, France (2016).[33] B. Schott, Stabilized cut finite element methods for complex interface coupled flow problems, Ph.D.thesis, Technical University of Munich (2017).[34] A. Massing, B. Schott, W. Wall, A stabilized Nitsche cut finite element method for the Oseenproblem, Comput. Methods Appl. Mech. Engrg. 328 (2018) 262–300.[35] M. Winter, B. Schott, A. Massing, W. Wall, A Nitsche cut finite element method for the Oseenproblem with general Navier boundary conditions, Comput. Methods Appl. Mech. Engrg. 330 (2017)220–252.[36] S. Zonca, L. Formaggia, C. Vergara, An unfitted formulation for the interaction of an incompressiblefluid with a thick structure via an xfem/dg approach, SIAM J. Sci. Comput. 40 (1) (2016) B59–B84.[37] L. Formaggia, C. Vergara, S. Zonca, Unfitted extended finite elements for composite grids, Comput.Math. Appl. 76 (4) (2018) 893–904.[38] T. E. Tezduyar, Computation of moving boundaries and interfaces and stabilization parameters,International Journal for Numerical Methods in Fluids 43 (5) (2003) 555–575.[39] J. Steger, F. Dougherty, J. Benek, A chimera grid scheme, in: Advances in Grid Generation, Vol.ASME FED-5, 1983, pp. 59–69.[40] J. L. Steger, J. A. Benek, On the use of composite grid schemes in computational aerodynamics,Computer Methods in Applied Mechanics and Engineering 64 (1-3) (1987) 301–320.[41] G. Chesshire, W. D. Henshaw, Composite overlapping meshes for the solution of partial differentialequations, J. Comput. Phys. 90 (1) (1990) 1–64.[42] G. Houzeaux, R. Codina, A Chimera method based on a Dirichlet/Neumann(Robin) coupling forthe Navier–Stokes equations, Comput. Methods Appl. Mech. Engrg. 192 (31-32) (2003) 3343–3377.[43] W. A. Wall, P. Gamnitzer, A. Gerstenberger, Fluid–structure interaction approaches on fixed gridsbased on two different domain decomposition ideas, International Journal of Computational FluidDynamics 22 (6) (2008) 411–427.[44] A. Hansbo, P. Hansbo, M. G. Larson, A finite element method on composite grids based on Nitsche’smethod, ESAIM: Mathematical Modelling and Numerical Analysis 37 (3) (2003) 495–514.[45] A. Massing, M. G. Larson, A. Logg, M. Rognes, A Nitsche-based cut finite element method for afluid-structure interaction problem, Commun. Appl. Math. Comput. Sci. 10 (2) (2015) 97–120.[46] B. Schott, S. Shahmiri, R. Kruse, W. Wall, A stabilized nitsche-type extended embedding mesh pproach for 3d low-and high-reynolds-number flows, Int. J. Numer. Methods Fluids.[47] A. Koblitz, S. Lovett, N. Nikiforakis, W. D. Henshaw, Direct numerical simulation of particulateflows with an overset grid method, Journal of Computational Physics 343 (2017) 414–431.[48] E. Burman, Projection stabilization of lagrange multipliers for the imposition of constraints oninterfaces and boundaries, Numerical Methods for Partial Differential Equations (2013) n/a–n/a.[49] E. Burman, P. Hansbo, Fictitious domain finite element methods using cut elements: I. A stabilizedLagrange multiplier method, Comput. Methods Appl. Mech. Engrg. 199 (2010) 2680–2686.[50] M. Puso, E. Kokko, R. Settgast, J. Sanders, B. Simpkins, B. Liu, An embedded mesh methodusing piecewise constant multipliers with stabilization: mathematical and numerical aspects, Int.J. Numer. Methods Eng.[51] A. Massing, M. G. Larson, A. Logg, M. E. Rognes, A stabilized Nitsche overlapping mesh methodfor the Stokes problem, Numerische Mathematik 128 (1) (2014) 73–101.[52] A. Johansson, M. G. Larson, A. Logg, High order cut finite element methods for the Stokes problem,Advanced Modeling and Simulation in Engineering Sciences 2 (1) (2015) 24.[53] J. Nitsche, ¨Uber ein Variationsprinzip zur L¨osung von Dirichlet-Problemen bei Verwendung vonTeilr¨aumen, die keinen Randbedingungen unterworfen sind. InAbhandlungen aus dem mathematis-chen Seminar der Universit¨at Hamburg 1971 Jul 1 (Vol. 36, No. 1, pp. 9-15).[54] J. M. Melenk, I. Babuˇska, The partition of unity finite element method: basic theory and applica-tions, Computer methods in applied mechanics and engineering 139 (1-4) (1996) 289–314.[55] Y. Huang, J. Xu, A conforming finite element method for overlapping and nonmatching grids, Math.Comp. 72 (243) (2003) 1057–1066.[56] A. Quarteroni, A. Valli, Numerical approximation of Partial Differential Equations, Springer, 1994,Ch. 9, pp. 297–337.[57] V. Girault, P.-A. Raviart, Finite element methods for Navier-Stokes equations: theory and algo-rithms, Vol. 5, Springer Science & Business Media, 2012.[58] C. Taylor, P. Hood, A numerical solution of the Navier-Stokes equations using the finite elementtechnique, Computers & Fluids 1 (1) (1973) 73–100.[59] A. Massing, M. G. Larson, A. Logg, Efficient implementation of finite element methods on non-matching and overlapping meshes in three dimensions, SIAM Journal on Scientific Computing 35 (1)(2013) C23–C47.[60] U. M. Mayer, A. Gerstenberger, W. A. Wall, Interface handling for three-dimensional higher-orderXFEM-computations in fluid–structure interaction, International Journal for Numerical Methodsin Engineering 79 (7) (2009) 846–869.[61] K. Wang, J. Gr´etarsson, A. Main, C. Farhat, Computational algorithms for tracking dynamic fluid–structure interfaces in embedded boundary methods, International Journal for Numerical Methodsin Fluids 70 (4) (2012) 515–535.[62] D. Nordsletten, N. Smith, D. Kay, A preconditioner for the finite element approximation to thearbitrary Lagrangian–Eulerian Navier–Stokes equations, SIAM Journal on Scientific Computing32 (2) (2010) 521–543.[63] A. Hessenthaler, O. R¨ohrle, D. Nordsletten, Validation of a non-conforming monolithic fluid-structure interaction method using phase-contrast MRI, International journal for numerical methodsin biomedical engineering 33 (8) (2017) e2845.[64] V. Shamanskii, A modification of Newton’s method, Ukrainian Mathematical Journal 19 (1) (1967)118–122.[65] J. H. Sp¨uhler, J. Jansson, N. Jansson, J. Hoffman, 3D Fluid-Structure Interaction Simulation ofAortic Valves Using a Unified Continuum ALE FEM Model, Frontiers in physiology 9.[66] D. Engwirda, Unstructured mesh methods for the Navier-Stokes equations, Undergraduate Thesis,School of Engineering, University of Sidney.[67] D. Engwirda, Locally optimal Delaunay-refinement and optimisation-based mesh generation, PhDthesis, School of Mathematics and Statistics, University of Sydney.[68] Sch¨afer, Michael and Turek, Stefan and Durst, Franz and Krause, Egon and Rannacher, Rolf, Bench-mark computations of laminar flow around a cylinder, in: Flow simulation with high-performancecomputers II, Springer, 1996, pp. 547–566.[69] M. Hadjicharalambous, J. Lee, N. P. Smith, D. A. Nordsletten, A displacement-based finite elementformulation for incompressible and nearly-incompressible cardiac mechanics, Computer methods inapplied mechanics and engineering 274 (2014) 213–236.[70] J. Lee, A. Cookson, I. Roy, E. Kerfoot, L. Asner, G. Vigueras, T. Sochi, S. Deparis, C. Michler,N. P. Smith, et al., Multiphysics computational modeling in CHeart, SIAM Journal on ScientificComputing 38 (3) (2016) C150–C178. igure 1: The problem domain in three instances: continuous, classical, single mesh setup and PUFEMsetup. In the later, the discrete domain is subdivided into the embedded and background overlappingmeshes. Both meshes contain a set of constrained nodes used to avoid ill-conditioning. Here we markthem with circles. Level 1 Level 2 Level 3 Level 4 Level 5Classich 0.1 0.05 0.025 0.0125 0 . / h . / h . / h . /
24 0 . / Table 1: Statistics corresponding to the mesh sets used in the steady Stokes and Navier-Stokes conver-gence tests.
FEM PUFEMLabel FEM1 PU1 PU2Background Embedded Background EmbeddedElements 14684 14816 1106 56362 4390DOF 67208 67745 5189 255742 20171 h (cm) 2.5 2.5 2.5 1.25 1.25 Table 2: Number of elements, DOF and average element size for the meshes used in the Turek benchmark. igure 2: (Left) An example of the embedded mesh with marked fluid-fluid and fluid-solid interfaces.(Right) The resulting weighting field ( ψ h )Figure 3: The intersection of one background element and two embedded elements at different stagesof computing the PUFEM weak form. (A) Identifying the intersection polygons. B Creating a newtessellation for the overlap area. (C)
Defining new set quadrature points for the sub-elements andre-evaluating the basis functions.
Time steps cd cl St ∆ p Lit. [68] - 3 . − .
24 0 . − .
01 0 . − .
305 2 . − . .
054 0 .
906 0 . . .
053 0 .
904 0 . . .
060 0 .
926 0 . . .
061 0 .
929 0 . . .
175 0 .
950 0 . . .
179 0 .
970 0 . . .
182 0 .
976 0 . . Table 3: Coefficient estimates from PUFEM and classic FEM for the Turek benchmark. igure 4: Domain used in the steady-state Stokes and Navier-Stokes problems.Figure 5: Example of solution fields for the Stokes problem obtained using the PUFEM and classicapproaches. The PUFEM figures include the background and embedded components, as well as thetotal (resulting from the weighted sum). Domain dimensions
Inflow width Max width Aortic length Sinus length Valve heightDistance (mm) 14 17.9 55.3 15.8 9.8
Fluid parameters Solid parameters µf ( Pa · s ) 3 × − µs ( Pa · s ) 20 × ρ ( kg/m
3) 1030
Table 4: (Top) The main domain dimensions in the aortic valve problem. (Bottom) The constitutivelaw parameters for the fluid and solid problems. igure 6: Plot showing the errors in the velocity and pressure fields for the Stokes flow problem obtainedwith the classical and PUFEM approaches. The dotted line indicates the optimal convergence rate.Figure 7: Example of solution fields for the Navier Stokes problem at Re
100 obtained using the PUFEMand classic approaches. The PUFEM figures include the background and embedded components, as wellas the total (resulting from the weighted sum). igure 8: The error plots for steady state incompressible Navier-Stokes. Both classic (C) and PUFEM(PU) methods are run at Reynolds number of 30 and 100.Figure 9: The domain of the Sch¨afer-Turek benchmark problem with boundary labels as well as examplemeshes used in PUFEM and classic approaches. igure 10: Compiled results for the Turek benchmark. (Top) A breakdown of the PUFEM solution atpeak c l (9.57 s in simulation time) into background, embedded and total fields for each component ofthe solution. The classic result is introduced for comparison. (Bottom) Velocity magnitude over a cyclein both PUFEM and classical simulations. Contour lines where added to aid comparison. igure 11: Compiled results for the ALE problem. (Top) A breakdown of the PUFEM solution intobackground(global), embedded and total. Classical result is introduced for comparison. (Middle) The x component of the velocity as it evolves between the two points of maximum displacement from the centre.Contour lines were added to aid comparison. (Bottom) Estimated drag force and pressure differenceacross the cylinder for the PUFEM and classical (BF1) methods. igure 12: An example of transient fixed nodes on the background mesh. The red arrows point to freenodes, while the blue ones point to their constrained counterpart at different time. The element colourcoding indicates the following: element cut by Γ h ff (deep grey), element completely weighed out by ψ h (light grey).Figure 13: (Left) Illustration of the fluid and valve domains used in the FSI test. A summary of thedomain dimensions is found in Table 4. (Right) The meshes used in the classic and PUFEM approaches.For visualization purposes the solid mesh was excluded. igure 14: (A) Six snapshots of the velocity magnitude as computed by the PUFEM and classicalapproaches. (B) The PUFEM and ALE mesh setups at the medium and minimum mesh deflection.At minimum deflection, the ALE mesh reaches a point of almost maximum deterioration. (C) Thedistribution in time of elements as a function of quality for both approaches. igure 15: The x and y displacements of the valve tip as it evolved in the classic and PUFEM approaches.The reference nodes where chosen such that they have identical coordinates at t = 0 s ..