A partition of unity approach to fluid mechanics and fluid-structure interaction

02/16/2019 ∙ by Maximilian Balmus, et al. ∙ NTNU KTH Royal Institute of Technology King's College London 0

For problems involving large deformations of thin structures, simulating fluid-structure interaction (FSI) remains challenging largely due to the need to balance computational feasibility, efficiency, and solution accuracy. Overlapping domain techniques have been introduced as a way to combine the fluid-solid mesh conformity, seen in moving-mesh methods, without the need for mesh smoothing or re-meshing, which is a core characteristic of fixed mesh approaches. In this work, we introduce a novel overlapping domain method based on a partition of unity approach. Unified function spaces are defined as a weighted sum of fields given on two overlapping meshes. The method is shown to achieve optimal convergence rates and to be stable for steady-state Stokes, Navier-Stokes, and ALE Navier-Stokes problems. Finally, we present results for FSI in the case of a 2D mock aortic valve simulation. These initial results point to the potential applicability of the method to a wide range of FSI applications, enabling boundary layer refinement and large deformations without the need for re-meshing or user-defined stabilization.



There are no comments yet.


This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Fluid-structure interaction (FSI) problems involving thin solids which undergo large deformations and translations can be encountered in a significant number of engineering applications. In industry, we have examples such as the design of parachutes stein1997parallel ; kalro2000parallel ; takizawa2011space and wind-turbines bazilevs20113d ; bazilevs20113d2 . In cardiovascular research, the simulation of valves hsu2015dynamic ; gao2017coupled and implanted devices mccormick2014computational holds a great potential for better understanding and treatment of a number of pathologies such as valve stenosis, regurgitation, heart failure and outflow obstruction. Building such models, however, remains challenging due to the need to balance computational costs and solution accuracy. In the case of cardiac valves, for example, studies typically require simplifications of the domain’s geometry lau2010mitral or the fluid models baillargeon2015human .

FSI approaches for this class of problems can be grouped into three main categories: interface-tracking, interface-capturing and overlapping domains methods bazilevs2013computational ; verkaik2014overlapping . In the case of interface-tracking, the fluid problem is typically based on the Arbitrary Lagrangian-Eulerian (ALE) hirt1974arbitrary ; donea1982arbitrary ; hughes1981lagrangian formulation. This allows for the fluid domain to deform with the solid and enables adjusting the element resolution close to surfaces in order to more accurately represent boundary layers. In van2007comparison , an ALE based method is shown to produce superior results when considering moderate deformations to three variations of fictitious domain (examples of interface-capturing) for similar mesh resolutions. However, for large deformations, it is known that the distortion of the fluid mesh can diminish the quality of elements and negatively impact the accuracy of solutions bazilevs2013computational . While multiple re-meshing techniques tezduyar2007modelling

have been proposed, the process introduces grid interpolation errors and its computational cost effectiveness is linked to the frequency at which the mesh needs to be adjusted. As shown in 

bavo2016fluid , where the particular case of cardiac valves is discussed, the ALE based simulations take more time to run than competing interface-capturing methods.

In contrast, interface-capturing methods do not require boundary fitted meshes for the fluid and solid and thus avoid the need for re-meshing. Examples include the Fictitious Domain Method glowinski1994fictitious ; glowinski1999distributed (FDM), where the coupling between the solid and fluid is achieved via additional Lagrange multiplier terms babuvska1973finite , and the Immersed Finite Element method (IFEM) zhang2004immersed ; ZhangGay2007 ; BoffiCavalliniGastaldi2015 ; BoffiGastaldi2017 , where the kinematic constraints between the fluid 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 Immersed Structural Potential Method (ISPM) gil2010immersed 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 result in poor approximations of pressure jumps and surface stresses dos2008partitioned ; kamensky2015immersogeometric . For this reason, mesh-adaptation van2004combined and XFEM enrichment gerstenberger2008extended ; alauzet2016nitsche ; Landajuela2016 ; Schott2017 ; MassingSchottWall2017 ; WinterSchottMassingEtAl2017 ; ZoncaFormaggiaVergara2016 ; FormaggiaVergaraZonca2018 techniques have been proposed in order to overcome these issues. Also a significant challenge for interface-capturing, as noted in tezduyar2003computation , remains the fact the resolution of the fluid flow around the structural surface is limited by the local element size of the fluid mesh. In practice, this leads to over-refinement of the fluid mesh along the moving trajectory of the solid and to significant increases of the computational costs.

Overlapping domain techniques StegerDoughertyBenek1983 ; StegerBenek1987 ; ChesshireHenshaw1990 ; HouzeauxCodina2003 ; WallGamnitzerGerstenberger2008 ; hansbo2003finite ; MassingLarsonLoggEtAl2015 ; SchottShahmiriKruseEtAl2016 ; verkaik2014overlapping ; KoblitzLovettNikiforakisEtAl2017 have been proposed with the aim of combining the advantages of interface-tracking and interface capturing techniques: fluid-solid mesh conformity, boundary layer tracking and eliminating re-meshing. The crux of these methods is the decomposition of the fluid problem into a background coarse component and solution-enriching embedded component that envelops the structures. A challenging aspect, however, is the coupling of the two fluid domains which has to be done weakly due the non-matching fluid-fluid interface. The use of Lagrange multipliers for example, see gerstenberger2008extended , is impeded by the need to properly choose function spaces in order to guarantee that the inf-sup condition holds for arbitrary moving interfaces. Alternatively, stabilization techniques can be employed to circumvent the inf-sup condition Burman2013 ; BurmanHansbo2010 ; PusoKokkoSettgastEtAl2014 . Overlapping mesh methods for the Stokes problem massing2014stabilized ; JohanssonLarsonLogg2015 ; SchottShahmiriKruseEtAl2016 which use Nitsche’s method nitschevariationsprinzip avoid introducing an additional Lagrange multiplier field, but nevertheless they require additional, parameter-dependent stabilization terms for the velocity and pressure jump in the vicinity of the fluid-fluid interface to guarantee optimal convergence rates and good system conditioning irrespective of the particular overlap configuration.

In this paper, we propose a new flexible and robust overlapping domain method which uses the partition of unity (PUFEM) approach melenk1996partition ; HuangXu2003 to decompose the fluid domain into a 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 using Taylor-Hood elements are defined. In the final finite element formulation of the fluid problem, a unified global function space is then used which is defined by taking a properly weighted sum of the function spaces associated with each mesh. To avoid ill-conditioning of the resulting system, additional constraints are introduced in (parts of) the overlap region.

The rest of this paper is structured as follows. After briefly recalling the classical mixed finite element approach for the Stokes problem in Section 2.1, we introduce its PUFEM based overlapping domain formulation including a detailed description of the 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 in Section 3, we explain how to combine the PUFEM approach with an ALE formulation of the Navier-Stokes problem to treat moving fluid domain and fluid–structure interaction problems. In Section 4, a number of numerical experiments are conducted. First, we investigate the stability and accuracy of our PU approach for a number of fluid flow problems posed on fixed static domains, see Section 4.14.3. To demonstrate the capability to handle large changes in the fluid domain geometry, we then consider a fluid flow driven by a oscillating cylinder in Section 4.4 before we turn to a full FSI problem in Section 4.5, where we compare a classical ALE-based approach with our novel combined ALE-PUFEM discretization for a two-dimensional mock aortic valve simulation. Finally in 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. We begin by reviewing the classic mixed FEM approach to Stokes flow (Section 2.1) and present 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 intersections of overlapping elements and perform the necessary integrations.

2.1 Classical mixed FEM approach to Stokes problems

Let be an arbitrary fluid domain on which we solve our problem. Since our focus is oriented towards FSI simulations, we also introduce a solid domain which for now, we consider to be fixed and rigid. The boundary of the fluid problem, , is composed of three, non-overlapping regions: the portions where Dirichlet and Neumann boundary conditions are applied, ( and , respectively), and the fluid-solid interface, . We then can write Stokes problem as: find the such that,


where is the fluid dynamic viscosity constant. For simplicity, we consider that the Neumann boundary condition tractions are null, i.e. .

In the case of the Stokes problem in (1), the classic continuous weak form can be written as: find such that,


where , and . Here the subscripts and indicate that the , subspaces are built such that they incorporate the Dirichlet and zero boundary value conditions on . Proofs of the well-posedness of this problem can be found in quarteroni1994introduction and girault2012finite . In the discrete setting, the resulting weak form is: find such that,


In this study, we will use the classic approach in (3) to compare with the PUFEM approach. In this case, we use the LBB stable Taylor-Hood elements taylor1973numerical . Thus, the discrete test and trial function spaces can be defined as:


Here is a discrete equivalent of , defined by a tessellation . denotes the element size defined as the diameter of the circumcircle. In preparation for the PUFEM discussion, we can also define as the discrete solid domain. designates the set of polynomial function of order .

Previous works quarteroni1994introduction ; girault2012finite derive a priori

error estimates where:


where and are positive constants independent of , and and denote the and norms, respectively. Using elements, if , then from interpolation theory quarteroni1994introduction :


2.2 PUFEM approach to Stokes problems

The core difference between the classic approach and the PUFEM setup is that the latter 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, we assume the background domain, , encompasses the entirety of the discrete fluid and solid domains. The embedded domain, , which satisfies , is designed to incorporate the solid and extends into fluid domains providing a boundary layer. The fluid and solid regions of are separated by the interface. Additionally, is the outer boundary of and serves as its interface with . The two meshes each have a corresponding tessellation, and , and element sizes denoted by and , respectively. Let and be the set of nodes defining the background mesh in the case of quadratic and linear interpolations. and are their embedded analogues. For the remainder of this paper, we consider the case where the meshes are 2D and triangular.

The discrete weak form for the PUFEM approach can be written as: find such that,


where , and are PUFEM function spaces that represent the counterparts of , and introduced in Section 2.1. The PUFEM function spaces are defined as the weighted sums of function spaces with support on the background and embedded meshes. Thus, based on the same notation of Melenk and Babuska melenk1996partition , the spaces can be written as:


where the local spaces are defined as:


for , , and


Note that with the exception of , all the discrete local spaces are continuous.

We allow the functions in to be discontinuous across in order to be able to represent the pressure jump between the fluid and embedded solid. The and sets contain nodes that are constrained in order to avoid ill-conditioning. More details on how these sets are constructed can be found in Section 2.2.1.

We define the field ( is the scalar equivalent of ) such that , its codomain is and on . For notational simplicity, we consider that all fields defined on the embedded mesh (i.e. , and ) are zero outside the region of overlap. In practice, we built using the following series of steps. First, we solved a diffusion problem on the embedded mesh, where we constrain the field to be zero on and ten on . Subsequently, we apply a upper boundary threshold such that the maximum node value is one. Finally, we smoothen the field using the hermitian polynomial . An example of the resulting field in the case of circular mesh can be seen Fig. 2. This approach was used irrespective of element size for consistency. Note however, that this process was chosen on an ad hoc basis. A full investigation 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 it guarantees a smooth transition across , irrespective of the mesh size used in either background or embedded mesh. As a result, this should always prevent having mismatched fluxes (i.e. when computed on both sides of the boundary) and the loss of mass across the interface.

2.2.1 Potential sources of ill-conditioning and constraints

Previously we introduced four sets of nodes that we want to constrain , , and and in this subsection we will discuss both their necessity and how we decide which nodes belong to them. With the introduction of the weighting field, the standard basis function support concept becomes insufficient to understand the contribution of each DOF to the total solution. Thus, we introduce a new metric of effective support fraction which for an embedded DOF we define as:



is the basis function of a degree of freedom corresponding to the

node coordinate. Note that in the case of a background DOF we replace with . As due to a particular overlap configuration, the contribution of DOFs corresponding to become very small. If then the contribution to the PUFEM solution is null and the system matrix becomes singular. Thus, we define the four sets of constrained nodes as follows:


In the case of the constraints on , i.e. (12a) and (12b), we have not observed particularly small values of the metric. However, we decided to include these nodes in order to obtain a smoother solution for the embedded mesh. refers to the area of the background domain comprised of elements completely covered by the embedded mesh. The definition in (12c) is meant to prevent any ill-conditioning, due to having a non-unique solution. The second condition, that the background nodes do not lie on fluid-fluid interface, is meant to prevent a circular definition (e.g. , therefore ). In the case of (12d), we chose to reduce the number of constrained DOF compared to and generally assigned to be 0.1. In tests which are not shown here, we saw that increasing the number of constrained pressure nodes can lead to a deterioration of the incompressibility condition at the domain level. While this effect can be lowered by decreasing , the trade off is an increase in the system’s condition number.

2.2.2 Fully discrete form of PUFEM for Stokes problems

Following the definition in Eq. (8), (10) and (9), the discrete fields and can be expanded into the sum of weighted basis functions:


and denote the piecewise polynomial functions defined on global and embedded meshes, respectively. Their subscripts are used in order to differentiate between first order () and second order () functions, while the superscript marks the degree of freedom (DOF) index. Thus, the four sets of basis functions can be written as: , , and . , , and are the total numbers of DOFs corresponding to the , , and

spaces. The nodal DOF vectors

and contain entries each (e.g. ). Note that, by design, this definition varies according to where we evaluate the field. Thus, for we have the definition in Equation (13a), while for the expansion simplifies to:

We can re-arrange the DOFs into four distinct vectors of unknown. For the background mesh, for example, we have:

The and vectors are defined analogously.

Algebraically, we can re-write the PUFEM weak form in Eq. (7) as:


Here the and blocks result from the Laplacian and incompressibility terms, respectively. Their general structure takes the form:

where and are placeholders for any and permutation. The inner sub-blocks are defined as:


Vectors and are the right hand side vectors resulting from imposed Dirichlet conditions. In order to apply the constraints on the system matrix, we need to distinguish between the Dirichlet condition and the constraints defined in Section 2.2.1. While the former is trivial, in the case of the latter, the original system matrix line is replaced with a set of local basis functions evaluated at the constrained coordinate. For example, if we wanted to constrain an arbitrary background velocity degree of freedom, , with the spatial coordinate , the new line of the system of equation would be equivalent to:


All other nodes, whether background or embedded, velocity or pressure, are fixed analogously.We use to denote the set of nodes of embedded element in which we find . We use to denote the set of nodes of embedded element in which we find . In the next section, we expound the process involve in computing the weighted weak form integral defined in  (15).

2.3 Computing the PUFEM weak form

We will start by illustrating this process with . The definition of this block, as defined in Eq. (15a) can be rewritten as the sum of a classic weak form matrix, defined on the background mesh, and a mixed classic-PUFEM matrix defined on the overlap:


Thus we can separate the construction of into two stages: first building a Laplacian block using classic FEM and second subtracting the classic weak form from the overlap area and replacing it with the PUFEM version. While the former is trivial, the later is more 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-elements which allows us relate the background and embedded meshes. Thus, at the start of the second stage, for each background element we want to identify all potentially intersecting embedded elements. This can be achieved optimally by adapting one of several approaches found in the literature massing2013efficient ; mayer2009interface ; wang2012computational . We then test the pairs of background and embedded candidates for intersection and we identify the area of overlap, see Fig. 3. If the area is non-zero, then the resulting intersection will be a convex polygon with up to 6 sides. Each area is subdivided into up to 4 subelements using a Delauney triangulation algorithm and these are stored in the tertiary mesh. The quality of this new mesh is not relevant since it is used for integration exclusively. To finalize the process, we loop over the intersection mesh and we update the block by subtracting the overlap classic weak form and replacing it with the PUFEM one.

Both and are defined exclusively on and thus do not require the first stage introduced . Since both are defined using fields from both meshes, the blocks are constructed by looping over the intersection mesh. In the case of , the block can built using the classic approach due to matching topologies for the space and fields.

The 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 application of PUFEM in the case of more complex settings, namely for ALE Navier-Stokes and FSI problems. Similarly to the Stokes flow case, we begin with a short review of the classic 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 a transient overlap between the two meshes. Finally, in Section 3.5 we outline a method of coupling the fluid solver with a hyper-elastic solid problem.

3.1 Classic FEM approach to the ALE Navier-Stokes problem

Let us consider a time dependent physical domain . At each time point , we use the notation to refer to the current spatial configuration of the fluid domain and let boundary be its boundary. We also define the reference domain which can be bijectively mapped to , for all . Hence, denotes a mapping function used to relate the reference and spatial domains. The non-conservative incompressible Navier-Stokes equation in arbitrary Lagrangian-Eulerian (ALE) form can be written as follows bazilevs2013computational ; nordsletten2010preconditioner : find such that


for all . Compared to the Stokes flow in (1

), the fluid model is augmented with a series of additional elements. The stress tensor is replaced by the symmetric Cauchy stress tensor:


Inertial effects are introduced using the time derivative and nonlinear advection terms. The operator is the time derivative with respect to a fixed point in . The ALE advection term is used to account for the arbitrary motion of the domain, where the arbitrary domain velocity field is defined as  nordsletten2010preconditioner ; bazilevs2013computational . Finally, we continue to assume that .

Using BDF(2) discretization, the classic non-conservative FEM problem for equation (18) at a given time step can be written as nordsletten2010preconditioner : find such that we have


Here, we consider the to be the velocity fields at times and , respectively, which have been transported via ALE mapping into the current spatial configuration.

Due to the non-linearity of the system, we approximate the solution using a Newton-Raphson (NR) based algorithm hessenthaler2017validation ; shamanskii1967modification . The problem is solved semi-monolithically in that the mesh velocity field, , is solved separately. If is provided analytically or is obtained from a weakly-coupled solid solver, then this process can be done at the beginning of each time step. If the fluid and solid solvers are coupled in a monolithic system, then the mesh velocity is recomputed prior to each NR iteration in order to account for the current estimate of the solid deformation.

3.2 Computing mesh velocity

The deformation of the mesh can be determined by an arbitrary problem which is guided by the known deformations on portions of the boundary. It is also well known that the choice of the problem can significantly impact mesh quality over time. For this reason, we propose to propagate the deformation of the boundary onto the entire mesh 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 we consider the material to be heterogeneous as we allow each element to stiffen in a squared inverse relationship with a metric of its relative element quality. A similar idea has been successfully implemented in spuhler20183d .

Let us consider the case where we want to compute the mesh deformation for time . We choose to be the reference undeformed domain. Thus the solid problem can be written as follows: find the arbitrary domain deformation field such that


This relates to the domain deformation velocity as . Here, the subscript indicates that the divergence operator is defined in the Lagrangian coordinates and the first Piola-Kirchhoff tensor is given as:


where is the deformation gradient and . The heterogeneous stiffness parameter , , is defined as:

Here is the reference stiffness; is an element quality metric defined as the ratio between the incircle () and circumcircle () of input triangular element; and denote the original and current shape of the element, respectively. As the quality of the element deteriorates, decreases, making it stiffer and, thus, less likely to further deform.

is a penalty factor on volumetric deformation. This choice was done 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 the reference, undeformed domain. For the normalization factor , we use the element at time zero as our reference. Given a discrete time point , we use the deformation field to compute the current space mapping function follows:


Based on this, the mesh velocity is defined as .

3.3 PUFEM approach to the ALE Navier-Stokes problem

In Section 2.2 we discussed the setup of PUFEM for flow around a rigid solid. Now we expand these original concepts in order to treat the case where the two meshes are allowed to move independently from each other.

A key requirement for the extension of PUFEM to ALE is the treatment of the material time derivative. Let and be the two reference domains for the background and embedded components. and are the two arbitrary mapping functions, where , are the respective physical domains. We also introduce the and as the partial time derivatives with respect to fixed points the global and embedded reference domains, respectively. For ease, we consider the ALE expansion of the material time derivative on an arbitrary semi-discrete (in space) scalar field which can be written in PUFEM form as , where , . Based on this, we can split the time material derivative of into three terms:

Furthermore, we choose to be the reference domain of and for and .

Thus, we can continue to expand the first term into its ALE form such that:

Similarly, we also get:

In the last case, we assumed that . The total material time derivative can thus be expressed as:

Moving to the discrete form, let us define and , the global and embedded meshes, respectively, for a given time step . Given that , , and change with the alignment of the meshes with respect to each other at each time step, let us define , and as the test and trial spaces at discrete time . Thus the new discrete weak form of the ALE Navier-Stokes problem can be written as follows: find such that for any we have that


where . Both partial derivatives and have been approximated using the BDF(2) scheme. In order to solve this system of equations, we again employ the technique described in shamanskii1967modification and hessenthaler2017validation . The structure of the resulting Jacobian matrix is largely based on that of system matrix described in Eq. (14). The changes arise in the sub-block where in addition to the Laplacian we also have mass and advection terms resulting from linearisation of the residual: